FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

octreeDataFaceList.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025 
00026 \*---------------------------------------------------------------------------*/
00027 
00028 #include "octreeDataFaceList.H"
00029 #include <meshTools/octree.H>
00030 
00031 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00032 
00033 defineTypeNameAndDebug(Foam::octreeDataFaceList, 0);
00034 
00035 Foam::scalar Foam::octreeDataFaceList::tol = 1E-6;
00036 
00037 
00038 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00039 
00040 void Foam::octreeDataFaceList::calcBb()
00041 {
00042     allBb_.setSize(faceLabels_.size());
00043     allBb_ = treeBoundBox
00044     (
00045         vector(GREAT, GREAT, GREAT),
00046         vector(-GREAT, -GREAT, -GREAT)
00047     );
00048 
00049     forAll (faceLabels_, faceLabelI)
00050     {
00051         label faceI = faceLabels_[faceLabelI];
00052 
00053         // Update bb of face
00054         treeBoundBox& myBb = allBb_[faceLabelI];
00055 
00056         const face& f = mesh_.localFaces()[faceI];
00057 
00058         forAll(f, fp)
00059         {
00060             const point& coord = mesh_.localPoints()[f[fp]];
00061 
00062             myBb.min() = min(myBb.min(), coord);
00063             myBb.max() = max(myBb.max(), coord);
00064         }
00065     }
00066 }
00067 
00068 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00069 
00070 // Construct from all faces in bMesh
00071 Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh)
00072 :
00073     mesh_(mesh),
00074     faceLabels_(mesh_.size()),
00075     allBb_(mesh_.size())
00076 {
00077     forAll(faceLabels_, faceI)
00078     {
00079         faceLabels_[faceI] = faceI;
00080     }
00081 
00082     // Generate tight fitting bounding box
00083     calcBb();
00084 }
00085 
00086 
00087 // Construct from selected faces in bMesh
00088 Foam::octreeDataFaceList::octreeDataFaceList
00089 (
00090     const bMesh& mesh,
00091     const labelList& faceLabels
00092 )
00093 :
00094     mesh_(mesh),
00095     faceLabels_(faceLabels),
00096     allBb_(faceLabels.size())
00097 {
00098     // Generate tight fitting bounding box
00099     calcBb();
00100 }
00101 
00102 
00103 
00104 // Construct as copy
00105 Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes)
00106 :
00107     mesh_(shapes.mesh()),
00108     faceLabels_(shapes.faceLabels()),
00109     allBb_(shapes.allBb())
00110 {}
00111 
00112 
00113 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
00114 
00115 Foam::octreeDataFaceList::~octreeDataFaceList()
00116 {}
00117 
00118 
00119 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00120 
00121 
00122 Foam::label Foam::octreeDataFaceList::getSampleType
00123 (
00124     const octree<octreeDataFaceList>& oc,
00125     const point& sample
00126 ) const
00127 {
00128     // Need to determine whether sample is 'inside' or 'outside'
00129     // Done by finding nearest face. This gives back a face which is
00130     // guaranteed to contain nearest point. This point can be
00131     // - in interior of face: compare to face normal
00132     // - on edge of face: compare to edge normal
00133     // - on point of face: compare to point normal
00134     // Unfortunately the octree does not give us back the intersection point
00135     // or where on the face it has hit so we have to recreate all that
00136     // information.
00137 
00138 
00139     // Find nearest face to sample
00140     treeBoundBox tightest(treeBoundBox::greatBox);
00141 
00142     scalar tightestDist = GREAT;
00143 
00144     label index = oc.findNearest(sample, tightest, tightestDist);
00145 
00146     if (index == -1)
00147     {
00148         FatalErrorIn
00149         (
00150             "octreeDataFaceList::getSampleType"
00151             "(octree<octreeDataFaceList>&, const point&)"
00152         )   << "Could not find " << sample << " in octree."
00153             << abort(FatalError);
00154     }
00155 
00156     label faceI = faceLabels_[index];
00157 
00158     // Get actual intersection point on face
00159 
00160     if (debug & 2)
00161     {
00162         Pout<< "getSampleType : sample:" << sample
00163             << " nearest face:" << faceI;
00164     }
00165 
00166     const face& f = mesh_.localFaces()[faceI];
00167 
00168     const pointField& points = mesh_.localPoints();
00169 
00170     pointHit curHit = f.nearestPoint(sample, points);
00171 
00172     //
00173     // 1] Check whether sample is above face
00174     //
00175 
00176     if (curHit.hit())
00177     {
00178         // Simple case. Compare to face normal.
00179 
00180         if (debug & 2)
00181         {
00182             Pout<< " -> face hit:" << curHit.hitPoint()
00183                 << " comparing to face normal " << mesh_.faceNormals()[faceI]
00184                 << endl;
00185         }
00186         return octree<octreeDataFaceList>::getVolType
00187         (
00188             mesh_.faceNormals()[faceI],
00189             sample - curHit.hitPoint()
00190         );
00191     }
00192 
00193     if (debug & 2)
00194     {
00195         Pout<< " -> face miss:" << curHit.missPoint();
00196     }
00197 
00198     //
00199     // 2] Check whether intersection is on one of the face vertices or
00200     //    face centre
00201     //
00202 
00203     // typical dimension as sqrt of face area.
00204     scalar typDim = sqrt(mag(f.normal(points))) + VSMALL;
00205 
00206     forAll(f, fp)
00207     {
00208         if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
00209         {
00210             // Face intersection point equals face vertex fp
00211 
00212             if (debug & 2)
00213             {
00214                     Pout<< " -> face point hit :" << points[f[fp]]
00215                         << " point normal:" << mesh_.pointNormals()[f[fp]]
00216                         << " distance:"
00217                         << mag(points[f[fp]] - curHit.missPoint())/typDim
00218                         << endl;
00219             }
00220             return octree<octreeDataFaceList>::getVolType
00221             (
00222                 mesh_.pointNormals()[f[fp]],
00223                 sample - curHit.missPoint()
00224             );
00225         }
00226     }
00227 
00228     // Get face centre
00229     point ctr(f.centre(points));
00230 
00231     if ((mag(ctr - curHit.missPoint())/typDim) < tol)
00232     {
00233         // Face intersection point equals face centre. Normal at face centre
00234         // is already average of face normals
00235 
00236         if (debug & 2)
00237         {
00238             Pout<< " -> centre hit:" << ctr
00239                 << " distance:"
00240                 << mag(ctr - curHit.missPoint())/typDim
00241                 << endl;
00242         }
00243 
00244         return octree<octreeDataFaceList>::getVolType
00245         (
00246             mesh_.faceNormals()[faceI],
00247             sample - curHit.missPoint()
00248         );
00249     }
00250 
00251 
00252     //
00253     // 3] Get the 'real' edge the face intersection is on
00254     //
00255 
00256     const labelList& myEdges = mesh_.faceEdges()[faceI];
00257 
00258     forAll(myEdges, myEdgeI)
00259     {
00260         const edge& e = mesh_.edges()[myEdges[myEdgeI]];
00261 
00262         pointHit edgeHit =
00263             line<point, const point&>
00264             (
00265                 points[e.start()],
00266                 points[e.end()]
00267             ).nearestDist(sample);
00268 
00269         point edgePoint;
00270         if (edgeHit.hit())
00271         {
00272             edgePoint = edgeHit.hitPoint();
00273         }
00274         else
00275         {
00276             edgePoint = edgeHit.missPoint();
00277         }
00278 
00279 
00280         if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
00281         {
00282             // Face intersection point lies on edge e
00283 
00284             // Calculate edge normal (wrong: uses face normals instead of
00285             // triangle normals)
00286             const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
00287 
00288             vector edgeNormal(vector::zero);
00289 
00290             forAll(myFaces, myFaceI)
00291             {
00292                 edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]];
00293             }
00294 
00295             if (debug & 2)
00296             {
00297                 Pout<< " -> real edge hit point:" << edgePoint
00298                     << " comparing to edge normal:" << edgeNormal
00299                     << endl;
00300             }
00301 
00302             // Found face intersection point on this edge. Compare to edge
00303             // normal
00304             return octree<octreeDataFaceList>::getVolType
00305             (
00306                 edgeNormal,
00307                 sample - curHit.missPoint()
00308             );
00309         }
00310     }
00311 
00312 
00313     //
00314     // 4] Get the internal edge (vertex - faceCentre) the face intersection
00315     //    is on
00316     //
00317 
00318     forAll(f, fp)
00319     {
00320         pointHit edgeHit =
00321             line<point, const point&>
00322             (
00323                 points[f[fp]],
00324                 ctr
00325             ).nearestDist(sample);
00326 
00327         point edgePoint;
00328         if (edgeHit.hit())
00329         {
00330             edgePoint = edgeHit.hitPoint();
00331         }
00332         else
00333         {
00334             edgePoint = edgeHit.missPoint();
00335         }
00336 
00337         if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
00338         {
00339             // Face intersection point lies on edge between two face triangles
00340 
00341             // Calculate edge normal as average of the two triangle normals
00342             label fpPrev = f.rcIndex(fp);
00343             label fpNext = f.fcIndex(fp);
00344 
00345             vector e = points[f[fp]] - ctr;
00346             vector ePrev = points[f[fpPrev]] - ctr;
00347             vector eNext = points[f[fpNext]] - ctr;
00348 
00349             vector nLeft = ePrev ^ e;
00350             nLeft /= mag(nLeft) + VSMALL;
00351 
00352             vector nRight = e ^ eNext;
00353             nRight /= mag(nRight) + VSMALL;
00354 
00355             if (debug & 2)
00356             {
00357                 Pout<< " -> internal edge hit point:" << edgePoint
00358                     << " comparing to edge normal "
00359                     << 0.5*(nLeft + nRight)
00360                     << endl;
00361             }
00362 
00363             // Found face intersection point on this edge. Compare to edge
00364             // normal
00365             return octree<octreeDataFaceList>::getVolType
00366             (
00367                 0.5*(nLeft + nRight),
00368                 sample - curHit.missPoint()
00369             );
00370         }
00371     }
00372 
00373     if (debug & 2)
00374     {
00375         Pout<< "Did not find sample " << sample
00376             << " anywhere related to nearest face " << faceI << endl
00377             << "Face:";
00378 
00379         forAll(f, fp)
00380         {
00381             Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
00382                 << endl;
00383         }
00384     }
00385 
00386     // Can't determine status of sample with respect to nearest face.
00387     // Either
00388     // - tolerances are wrong. (if e.g. face has zero area)
00389     // - or (more likely) surface is not closed.
00390 
00391     return octree<octreeDataFaceList>::UNKNOWN;
00392 }
00393 
00394 
00395 bool Foam::octreeDataFaceList::overlaps
00396 (
00397     const label index,
00398     const treeBoundBox& sampleBb
00399 ) const
00400 {
00401     return sampleBb.overlaps(allBb_[index]);
00402 }
00403 
00404 
00405 bool Foam::octreeDataFaceList::contains
00406 (
00407     const label,
00408     const point&
00409 ) const
00410 {
00411     notImplemented
00412     (
00413         "octreeDataFaceList::contains(const label, const point&)"
00414     );
00415     return false;
00416 }
00417 
00418 
00419 bool Foam::octreeDataFaceList::intersects
00420 (
00421     const label index,
00422     const point& start,
00423     const point& end,
00424     point& intersectionPoint
00425 ) const
00426 {
00427     label faceI = faceLabels_[index];
00428 
00429     const face& f = mesh_.localFaces()[faceI];
00430 
00431     const vector dir(end - start);
00432 
00433     // Disable picking up intersections behind us.
00434     scalar oldTol = intersection::setPlanarTol(0.0);
00435 
00436     pointHit inter =
00437         f.ray
00438         (
00439             start,
00440             dir,
00441             mesh_.localPoints(),
00442             intersection::HALF_RAY,
00443             intersection::VECTOR
00444         );
00445 
00446     intersection::setPlanarTol(oldTol);
00447 
00448     if (inter.hit() && inter.distance() <= mag(dir))
00449     {
00450         intersectionPoint = inter.hitPoint();
00451 
00452         return true;
00453     }
00454     else
00455     {
00456         return false;
00457     }
00458 }
00459 
00460 
00461 bool Foam::octreeDataFaceList::findTightest
00462 (
00463     const label index,
00464     const point& sample,
00465     treeBoundBox& tightest
00466 ) const
00467 {
00468     // Get nearest and furthest away vertex
00469     point myNear, myFar;
00470     allBb_[index].calcExtremities(sample, myNear, myFar);
00471 
00472     const point dist = myFar - sample;
00473     scalar myFarDist = mag(dist);
00474 
00475     point tightestNear, tightestFar;
00476     tightest.calcExtremities(sample, tightestNear, tightestFar);
00477 
00478     scalar tightestFarDist = mag(tightestFar - sample);
00479 
00480     if (tightestFarDist < myFarDist)
00481     {
00482         // Keep current tightest.
00483         return false;
00484     }
00485     else
00486     {
00487         // Construct bb around sample and myFar
00488         const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
00489 
00490         tightest.min() = sample - dist2;
00491         tightest.max() = sample + dist2;
00492 
00493         return true;
00494     }
00495 }
00496 
00497 
00498 // Determine numerical value of sign of sample compared to shape at index
00499 Foam::scalar Foam::octreeDataFaceList::calcSign
00500 (
00501     const label index,
00502     const point& sample,
00503     vector&
00504 ) const
00505 {
00506     label faceI = faceLabels_[index];
00507 
00508     const face& f = mesh_.localFaces()[faceI];
00509 
00510     point ctr = f.centre(mesh_.localPoints());
00511 
00512     vector vec = sample - ctr;
00513 
00514     vec /= mag(vec) + VSMALL;
00515 
00516     return mesh_.faceNormals()[faceI] & vec;
00517 }
00518 
00519 
00520 // Calculate nearest point on/in shapei
00521 Foam::scalar Foam::octreeDataFaceList::calcNearest
00522 (
00523     const label index,
00524     const point& sample,
00525     point& nearest
00526 ) const
00527 {
00528     label faceI = faceLabels_[index];
00529 
00530     const face& f = mesh_.localFaces()[faceI];
00531 
00532     pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints());
00533 
00534     if (nearHit.hit())
00535     {
00536         nearest = nearHit.hitPoint();
00537     }
00538     else
00539     {
00540         nearest = nearHit.missPoint();
00541     }
00542 
00543     if (debug & 1)
00544     {
00545         point ctr = f.centre(mesh_.localPoints());
00546 
00547         scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest);
00548 
00549         Pout<< "octreeDataFaceList::calcNearest : "
00550             << "sample:" << sample
00551             << "  faceI:" << faceI
00552             << "  ctr:" << ctr
00553             << "  sign:" << sign
00554             << "  nearest point:" << nearest
00555             << "  distance to face:" << nearHit.distance()
00556             << endl;
00557     }
00558     return nearHit.distance();
00559 }
00560 
00561 
00562 void Foam::octreeDataFaceList::write
00563 (
00564     Ostream& os,
00565     const label index
00566 ) const
00567 {
00568     os << faceLabels_[index] << " " << allBb_[index];
00569 }
00570 
00571 
00572 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines