BUGFIX: Uninitialised member variables
[foam-extend-3.2.git] / applications / utilities / surface / surfaceToPatch / treeDataPrimitivePatch.C
blob16d79e85f3faa9110207b8706f9e077a6a0fe4a1
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright held by original author
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
11     OpenFOAM is free software; you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by the
13     Free Software Foundation; either version 2 of the License, or (at your
14     option) any later version.
16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19     for more details.
21     You should have received a copy of the GNU General Public License
22     along with OpenFOAM; if not, write to the Free Software Foundation,
23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "treeDataPrimitivePatch.H"
28 #include "indexedOctree.H"
29 #include "triangleFuncs.H"
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 template
35     class Face,
36     template<class> class FaceList,
37     class PointField,
38     class PointType
40 Foam::scalar
41 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
42 tolSqr = sqr(1E-6);
45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
47 template
49     class Face,
50     template<class> class FaceList,
51     class PointField,
52     class PointType
54 Foam::treeBoundBox
55 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
56 calcBb
58     const pointField& points,
59     const face& f
62     treeBoundBox bb(points[f[0]], points[f[0]]);
64     for (label fp = 1; fp < f.size(); fp++)
65     {
66         const point& p = points[f[fp]];
68         bb.min() = min(bb.min(), p);
69         bb.max() = max(bb.max(), p);
70     }
71     return bb;
75 template
77     class Face,
78     template<class> class FaceList,
79     class PointField,
80     class PointType
82 void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
83 update()
85     if (cacheBb_)
86     {
87         bbs_.setSize(patch_.size());
89         forAll(patch_, i)
90         {
91             bbs_[i] = calcBb(patch_.points(), patch_[i]);
92         }
93     }
97 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
99 // Construct from components
100 template
102     class Face,
103     template<class> class FaceList,
104     class PointField,
105     class PointType
107 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
108 treeDataPrimitivePatch
110     const bool cacheBb,
111     const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
114     patch_(patch),
115     cacheBb_(cacheBb)
117     update();
121 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
123 template
125     class Face,
126     template<class> class FaceList,
127     class PointField,
128     class PointType
130 Foam::pointField
131 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
132 points() const
134     pointField cc(patch_.size());
136     forAll(patch_, i)
137     {
138         cc[i] = patch_[i].centre(patch_.points());
139     }
141     return cc;
145 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
146 //  Only makes sense for closed surfaces.
147 template
149     class Face,
150     template<class> class FaceList,
151     class PointField,
152     class PointType
154 Foam::label
155 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
156 getVolumeType
158     const indexedOctree
159     <
160         treeDataPrimitivePatch
161         <
162             Face,
163             FaceList,
164             PointField,
165             PointType
166         >
167     >& oc,
168     const point& sample
169 ) const
171     // Need to determine whether sample is 'inside' or 'outside'
172     // Done by finding nearest face. This gives back a face which is
173     // guaranteed to contain nearest point. This point can be
174     // - in interior of face: compare to face normal
175     // - on edge of face: compare to edge normal
176     // - on point of face: compare to point normal
177     // Unfortunately the octree does not give us back the intersection point
178     // or where on the face it has hit so we have to recreate all that
179     // information.
182     // Find nearest face to sample
183     pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
185     if (info.index() == -1)
186     {
187         FatalErrorIn
188         (
189             "treeDataPrimitivePatch::getSampleType"
190             "(indexedOctree<treeDataPrimitivePatch>&, const point&)"
191         )   << "Could not find " << sample << " in octree."
192             << abort(FatalError);
193     }
196     // Get actual intersection point on face
197     label faceI = info.index();
199     if (debug & 2)
200     {
201         Pout<< "getSampleType : sample:" << sample
202             << " nearest face:" << faceI;
203     }
205     const pointField& points = patch_.localPoints();
206     const face& f = patch_.localFaces()[faceI];
208     // Retest to classify where on face info is. Note: could be improved. We
209     // already have point.
211     pointHit curHit = f.nearestPoint(sample, points);
212     const vector area = f.normal(points);
213     const point& curPt = curHit.rawPoint();
215     //
216     // 1] Check whether sample is above face
217     //
219     if (curHit.hit())
220     {
221         // Nearest point inside face. Compare to face normal.
223         if (debug & 2)
224         {
225             Pout<< " -> face hit:" << curPt
226                 << " comparing to face normal " << area << endl;
227         }
228         return indexedOctree<treeDataPrimitivePatch>::getSide
229         (
230             area,
231             sample - curPt
232         );
233     }
235     if (debug & 2)
236     {
237         Pout<< " -> face miss:" << curPt;
238     }
240     //
241     // 2] Check whether intersection is on one of the face vertices or
242     //    face centre
243     //
245     const scalar typDimSqr = mag(area) + VSMALL;
247     forAll(f, fp)
248     {
249         if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
250         {
251             // Face intersection point equals face vertex fp
253             // Calculate point normal (wrong: uses face normals instead of
254             // triangle normals)
256             return indexedOctree<treeDataPrimitivePatch>::getSide
257             (
258                 patch_.pointNormals()[f[fp]],
259                 sample - curPt
260             );
261         }
262     }
264     const point fc(f.centre(points));
266     if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
267     {
268         // Face intersection point equals face centre. Normal at face centre
269         // is already average of face normals
271         if (debug & 2)
272         {
273             Pout<< " -> centre hit:" << fc
274                 << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
275         }
277         return indexedOctree<treeDataPrimitivePatch>::getSide
278         (
279             area,
280             sample - curPt
281         );
282     }
286     //
287     // 3] Get the 'real' edge the face intersection is on
288     //
290     const labelList& fEdges = patch_.faceEdges()[faceI];
292     forAll(fEdges, fEdgeI)
293     {
294         label edgeI = fEdges[fEdgeI];
295         const edge& e = patch_.edges()[edgeI];
297         pointHit edgeHit = e.line(points).nearestDist(sample);
299         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
300         {
301             // Face intersection point lies on edge e
303             // Calculate edge normal (wrong: uses face normals instead of
304             // triangle normals)
305             const labelList& eFaces = patch_.edgeFaces()[edgeI];
307             vector edgeNormal(vector::zero);
309             forAll(eFaces, i)
310             {
311                 edgeNormal += patch_.faceNormal()[eFaces[i]];
312             }
314             if (debug & 2)
315             {
316                 Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
317                     << " comparing to edge normal:" << edgeNormal
318                     << endl;
319             }
321             // Found face intersection point on this edge. Compare to edge
322             // normal
323             return indexedOctree<treeDataPrimitivePatch>::getSide
324             (
325                 edgeNormal,
326                 sample - curPt
327             );
328         }
329     }
332     //
333     // 4] Get the internal edge the face intersection is on
334     //
336     forAll(f, fp)
337     {
338         pointHit edgeHit = linePointRef
339         (
340             points[f[fp]],
341             fc
342         ).nearestDist(sample);
344         if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
345         {
346             // Face intersection point lies on edge between two face triangles
348             // Calculate edge normal as average of the two triangle normals
349             vector e = points[f[fp]] - fc;
350             vector ePrev = points[f[f.rcIndex(fp)]] - fc;
351             vector eNext = points[f[f.fcIndex(fp)]] - fc;
353             vector nLeft = ePrev ^ e;
354             nLeft /= mag(nLeft) + VSMALL;
356             vector nRight = e ^ eNext;
357             nRight /= mag(nRight) + VSMALL;
359             if (debug & 2)
360             {
361                 Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
362                     << " comparing to edge normal "
363                     << 0.5*(nLeft + nRight)
364                     << endl;
365             }
367             // Found face intersection point on this edge. Compare to edge
368             // normal
369             return indexedOctree<treeDataPrimitivePatch>::getSide
370             (
371                 0.5*(nLeft + nRight),
372                 sample - curPt
373             );
374         }
375     }
377     if (debug & 2)
378     {
379         Pout<< "Did not find sample " << sample
380             << " anywhere related to nearest face " << faceI << endl
381             << "Face:";
383         forAll(f, fp)
384         {
385             Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
386                 << endl;
387         }
388     }
390     // Can't determine status of sample with respect to nearest face.
391     // Either
392     // - tolerances are wrong. (if e.g. face has zero area)
393     // - or (more likely) surface is not closed.
395     return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
399 // Check if any point on shape is inside cubeBb.
400 template
402     class Face,
403     template<class> class FaceList,
404     class PointField,
405     class PointType
407 bool
408 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
409 overlaps
411     const label index,
412     const treeBoundBox& cubeBb
413 ) const
415     // 1. Quick rejection: bb does not intersect face bb at all
416     if (cacheBb_)
417     {
418         if (!cubeBb.overlaps(bbs_[index]))
419         {
420             return false;
421         }
422     }
423     else
424     {
425         if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
426         {
427             return false;
428         }
429     }
432     // 2. Check if one or more face points inside
434     const pointField& points = patch_.points();
435     const face& f = patch_[index];
437     forAll(f, fp)
438     {
439         if (cubeBb.contains(points[f[fp]]))
440         {
441             return true;
442         }
443     }
445     // 3. Difficult case: all points are outside but connecting edges might
446     // go through cube. Use triangle-bounding box intersection.
447     const point fc = f.centre(points);
449     forAll(f, fp)
450     {
451         bool triIntersects = triangleFuncs::intersectBb
452         (
453             points[f[fp]],
454             points[f[f.fcIndex(fp)]],
455             fc,
456             cubeBb
457         );
459         if (triIntersects)
460         {
461             return true;
462         }
463     }
464     return false;
468 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
469 // nearestPoint.
470 template
472     class Face,
473     template<class> class FaceList,
474     class PointField,
475     class PointType
477 void
478 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
479 findNearest
481     const labelList& indices,
482     const point& sample,
484     scalar& nearestDistSqr,
485     label& minIndex,
486     point& nearestPoint
487 ) const
489     const pointField& points = patch_.points();
491     forAll(indices, i)
492     {
493         label index = indices[i];
495         const face& f = patch_[index];
497         pointHit nearHit = f.nearestPoint(sample, points);
498         scalar distSqr = sqr(nearHit.distance());
500         if (distSqr < nearestDistSqr)
501         {
502             nearestDistSqr = distSqr;
503             minIndex = index;
504             nearestPoint = nearHit.rawPoint();
505         }
506     }
510 template
512     class Face,
513     template<class> class FaceList,
514     class PointField,
515     class PointType
517 bool
518 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
519 intersects
521     const label index,
522     const point& start,
523     const point& end,
524     point& intersectionPoint
525 ) const
527     // Do quick rejection test
528     if (cacheBb_)
529     {
530         const treeBoundBox& faceBb = bbs_[index];
532         if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
533         {
534             // start and end in same block outside of faceBb.
535             return false;
536         }
537     }
539     const pointField& points = patch_.points();
540     const face& f = patch_[index];
541     const point fc = f.centre(points);
542     const vector dir(end - start);
544     pointHit inter = patch_[index].intersection
545     (
546         start,
547         dir,
548         fc,
549         points,
550         intersection::HALF_RAY
551     );
553     if (inter.hit() && inter.distance() <= 1)
554     {
555         // Note: no extra test on whether intersection is in front of us
556         // since using half_ray
557         intersectionPoint = inter.hitPoint();
558         return true;
559     }
560     else
561     {
562         return false;
563     }
567 // ************************************************************************* //