1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright held by original author
7 -------------------------------------------------------------------------------
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
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 * * * * * * * * * * * * * //
36 template<class> class FaceList,
41 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
50 template<class> class FaceList,
55 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
58 const pointField& points,
62 treeBoundBox bb(points[f[0]], points[f[0]]);
64 for (label fp = 1; fp < f.size(); fp++)
66 const point& p = points[f[fp]];
68 bb.min() = min(bb.min(), p);
69 bb.max() = max(bb.max(), p);
78 template<class> class FaceList,
82 void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
87 bbs_.setSize(patch_.size());
91 bbs_[i] = calcBb(patch_.points(), patch_[i]);
97 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
99 // Construct from components
103 template<class> class FaceList,
107 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
108 treeDataPrimitivePatch
111 const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
121 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
126 template<class> class FaceList,
131 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
134 pointField cc(patch_.size());
138 cc[i] = patch_[i].centre(patch_.points());
145 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
146 // Only makes sense for closed surfaces.
150 template<class> class FaceList,
155 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
160 treeDataPrimitivePatch
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
182 // Find nearest face to sample
183 pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
185 if (info.index() == -1)
189 "treeDataPrimitivePatch::getSampleType"
190 "(indexedOctree<treeDataPrimitivePatch>&, const point&)"
191 ) << "Could not find " << sample << " in octree."
192 << abort(FatalError);
196 // Get actual intersection point on face
197 label faceI = info.index();
201 Pout<< "getSampleType : sample:" << sample
202 << " nearest face:" << faceI;
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();
216 // 1] Check whether sample is above face
221 // Nearest point inside face. Compare to face normal.
225 Pout<< " -> face hit:" << curPt
226 << " comparing to face normal " << area << endl;
228 return indexedOctree<treeDataPrimitivePatch>::getSide
237 Pout<< " -> face miss:" << curPt;
241 // 2] Check whether intersection is on one of the face vertices or
245 const scalar typDimSqr = mag(area) + VSMALL;
249 if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
251 // Face intersection point equals face vertex fp
253 // Calculate point normal (wrong: uses face normals instead of
256 return indexedOctree<treeDataPrimitivePatch>::getSide
258 patch_.pointNormals()[f[fp]],
264 const point fc(f.centre(points));
266 if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
268 // Face intersection point equals face centre. Normal at face centre
269 // is already average of face normals
273 Pout<< " -> centre hit:" << fc
274 << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
277 return indexedOctree<treeDataPrimitivePatch>::getSide
287 // 3] Get the 'real' edge the face intersection is on
290 const labelList& fEdges = patch_.faceEdges()[faceI];
292 forAll(fEdges, fEdgeI)
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)
301 // Face intersection point lies on edge e
303 // Calculate edge normal (wrong: uses face normals instead of
305 const labelList& eFaces = patch_.edgeFaces()[edgeI];
307 vector edgeNormal(vector::zero);
311 edgeNormal += patch_.faceNormal()[eFaces[i]];
316 Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
317 << " comparing to edge normal:" << edgeNormal
321 // Found face intersection point on this edge. Compare to edge
323 return indexedOctree<treeDataPrimitivePatch>::getSide
333 // 4] Get the internal edge the face intersection is on
338 pointHit edgeHit = linePointRef
342 ).nearestDist(sample);
344 if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
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;
361 Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
362 << " comparing to edge normal "
363 << 0.5*(nLeft + nRight)
367 // Found face intersection point on this edge. Compare to edge
369 return indexedOctree<treeDataPrimitivePatch>::getSide
371 0.5*(nLeft + nRight),
379 Pout<< "Did not find sample " << sample
380 << " anywhere related to nearest face " << faceI << endl
385 Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
390 // Can't determine status of sample with respect to nearest face.
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.
403 template<class> class FaceList,
408 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
412 const treeBoundBox& cubeBb
415 // 1. Quick rejection: bb does not intersect face bb at all
418 if (!cubeBb.overlaps(bbs_[index]))
425 if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
432 // 2. Check if one or more face points inside
434 const pointField& points = patch_.points();
435 const face& f = patch_[index];
439 if (cubeBb.contains(points[f[fp]]))
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);
451 bool triIntersects = triangleFuncs::intersectBb
454 points[f[f.fcIndex(fp)]],
468 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
473 template<class> class FaceList,
478 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
481 const labelList& indices,
484 scalar& nearestDistSqr,
489 const pointField& points = patch_.points();
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)
502 nearestDistSqr = distSqr;
504 nearestPoint = nearHit.rawPoint();
513 template<class> class FaceList,
518 Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
524 point& intersectionPoint
527 // Do quick rejection test
530 const treeBoundBox& faceBb = bbs_[index];
532 if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
534 // start and end in same block outside of faceBb.
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
550 intersection::HALF_RAY
553 if (inter.hit() && inter.distance() <= 1)
555 // Note: no extra test on whether intersection is in front of us
556 // since using half_ray
557 intersectionPoint = inter.hitPoint();
567 // ************************************************************************* //