[MMC15]

Hi OpenFOAMers!!

I'm looking for an algorithm to search the intersection between a line and cell faces of a mesh. I know Ray Tracing, which is quite similar to my problem, is an important field of study for graphic programmers, but I didn't find nothing practically interesting about CFD mesh applications. I've already searched in the forum but I didn't find nothing.

Could anybody help me?
Thank you a lot!

Regards from Italy!

[mirx]

There are a lot of functions in OF to find intersections of faces and lines/rays. Usually they are methods of indexedOctree classes. It only depends on, what you would like to intersect.
Boundary faces, triSurface faces, mesh faces?
With lines (find intersection between start point and end point) or with rays (start point and vector)?

Do you just need to find these intersections or are you really looking for a ray tracing algorithm, like monte carlo ray tracing?

[MMC15]

I'd like to know the intersection points of a line from a start point to an end point with every mesh faces at runtime.
The indexedOctree functions may be useful, but I'm having troubles with their use.
Have you some advices about their application?
Which headers have I to include?
Is there a specific function that you suggest?

[mirx]

Here are some code snippets I used to find all intersections of a line with all mesh faces. This is no complete working code.

#include "polyMesh.H"
#include "meshTools.H"
#include "treeDataFace.H"

    treeBoundBox allBb(mesh.points());


    scalar bbTol = 1e-6 * allBb.avgDim();


    point& bbMin = allBb.min();
    bbMin.x() -= bbTol;
    bbMin.y() -= bbTol;
    bbMin.z() -= bbTol;


    point& bbMax = allBb.max();
    bbMax.x() += 2*bbTol;
    bbMax.y() += 2*bbTol;
    bbMax.z() += 2*bbTol;


    indexedOctree<treeDataFace> faceTree
    (
        treeDataFace(false, mesh_),
        allBb, // overall search domain
        8, // maxLevel
        10, // leafsize
        3.0 // duplicity
    );

        const point pStart = ...; // start point of line

        const point pEnd = ...; // end point of line

        const vector eVec(pEnd - pStart); // line vector

        const scalar eMag = mag(eVec); // edge length
        const vector tolVec = 1e-6*eVec;        
        point p0 = pStart - tolVec;
        const point p1 = pEnd + tolVec;
        


        while(true)
        {
            pointIndexHit pHit = faceTree.findLine(p0, p1);
            
            if (pHit.hit())
            {
                const label faceI = pHit.index(); // face label of hit face

                const point hitPoint = pHit.hitPoint(); // intersection point

                

                const vector& area = mesh_.faceAreas()[pHit.index()];                scalar typDim = Foam::sqrt(mag(area));
                
                // stop search if new start point is near to edge end

                if ((mag(pHit.hitPoint() - pEnd)/typDim) < SMALL)
                {
                    break;
                }

                // set new start point shortly after previous start point

                p0 = pHit.hitPoint() + tolVec;
            }
            else
            {
                // No hit.
                break;
            }
        }

EDIT: The while loop is used to find ALL intersections. The faceTree.findLine(p0, p1) method only returns the first intersection starting from p0. To find all intersections you have to change the starting point, if findLine() found an intersection.

[MMC15]

Thank you mirx.
Nice code, it works very well, I just add a list to save the hit points.

DynamicList<pointIndexHit> pHitList;

pHitList.append(pHit);

I read that exists an other indexedOctree function called findLineAny, is there a reason why don't you use it, instead of findLine?

[mirx]

There was a reason that I didn't use it. I wanted to run some extra code for every intersection I had found.

But you are right. findLineAny actually executes a similar while loop as the one I have posted. I think my loop was inspired by the one in findLineAny. So if you just need the intersections, I think you can use something like:

const point pStart = ...; 
const point pEnd = ...; 
List<pointIndexHit> pHit(faceTree.findLineAny(pStart, pEnd));

[MMC15]

I tried that but the list output seems to be not accepted:

 error: no matching function for call to ‘Foam::List<Foam::PointIndexHit<Foam::Vector<double> > >::List(Foam::pointIndexHit)’
     List<pointIndexHit> pHit(faceTree.findLineAny(pStart, pEnd));’

[mirx]

Oh sorry, my mistake. I confused Any with All. triSurfaceSearch has a method called fineLineAll.

void findLineAll        
(
     const point& start,
     const point& end,
     List<pointIndexHit>&
) const;

This returns a list of all pointIndexHit. indexedOctree doesn't have this method. (Of course... triSurfaceSearch.findLineAll() uses indexedOctree's findLine function with a loop)

findLineAny returns any of the intersections, not all.

[MMC15]

OK! Thank you a lot for your useful replies!

Bye!

P.S.: I said "Any", so it was my mistake! :-)

[Yahoo]

Hey guys,

Is any body still interested in this ray tracing application in OpenFOAM?

I am having trouble with including the header file "meshtools.H"--fatal error, no such file or directory.

Do you have any comments?

[alexeym]

Hi,

it's meshTools.H not meshtools.h. Also you need

-I$(LIB_SRC)/meshTools/lnInclude

in your options file.

[Yahoo]

Thanks Alexeym!
I got it going!
one more question: do you know how can i use findLineAll member function of triSurfaceSearch?

[alexeym]

Hi,

Looked through the documentation on the method, now I have doubts I've got your question right. You've got triSurface, you construct triSurfaceSearch, you call findLineAll. So what's your question?

[mirx]

For a ray tracing application I think you will be interested in the first intersection of a ray with a surface. In this case you may use method findLine of indexedOctree<treeDataTriSurface>.

This is a short example how to use this method:

const triSurface surf(runTime.constantPath()/"triSurface"/surfFileName);
const vectorField& normals = surf.faceNormals();
const triSurfaceSearch querySurf(surf);
const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
point searchStart(0, 0, 0);
point searchEnd(1, 0, 0);
pointIndexHit pHit = tree.findLine(searchStart, searchEnd);
if ( pHit.hit() )
{
    point hitPoint = pHit.hitPoint();
    label hitTriangle = pHit.index();
    vector normal = normals[hitTriangle];
}

With the triangle label you are able to identify the plane normal and thus the interaction of the ray with the surface (reflection, refraction).

[sinsni]

Greetings! Hope someone will be able to give some hints on my issue. Thanks in advance.

I'm doing ray tracing, using octree:

indexedOctree<treeDataFace> faceTree     
(         
     treeDataFace(false, mesh),         
     allBb, // overall search domain         
     8, // maxLevel         
     10, // leafsize         
     3.0 // duplicity     
);

and then

faceTree.findLine(p0, p1);

I have a mesh with millions of cells. With mesh this big, performance of the octree search becomes quite poor. But I don't need to build indexedOctree from the whole mesh, like I do:

treeDataFace(false, mesh),

All my rays are in a restricted part of the mesh. Is it possible to do kind of like this:

treeDataFace(false, "small part of the mesh"),

Ultimately, I want to improve performance. Maybe there's some other way...

[David*]

Originally Posted by mirx

For a ray tracing application I think you will be interested in the first intersection of a ray with a surface. In this case you may use method findLine of indexedOctree<treeDataTriSurface>.

This is a short example how to use this method:

const triSurface surf(runTime.constantPath()/"triSurface"/surfFileName);
const vectorField& normals = surf.faceNormals();
const triSurfaceSearch querySurf(surf);
const indexedOctree<treeDataTriSurface>& tree = querySurf.tree();
point searchStart(0, 0, 0);
point searchEnd(1, 0, 0);
pointIndexHit pHit = tree.findLine(searchStart, searchEnd);
if ( pHit.hit() )
{
    point hitPoint = pHit.hitPoint();
    label hitTriangle = pHit.index();
    vector normal = normals[hitTriangle];
}

...

Small note to save others the amounts of time I spent on this: To use this (very useful!) code without the wmake script (i.e., your own Makefile), it is necessary to apply the -DNoRepository compiler flag. Otherwise, you'll run into something like

undefined reference to `Foam::PrimitivePatch<Foam::labelledTri, Foam::List, Foam::Field<Foam::Vector<double> >, Foam::Vector<double> >::faceNormals() const'

because PrimitivePatch.H is a template class and does not automatically include the necessary .C files.

Cheers, David

[tom_flint2012]

Sinsi,

Did you ever find a method of restricting the search to a smaller part of the mesh?

All the best,

Tom
treeDataFace(false, "small part of the mesh"),

[ashish.svm]

is it possible to search the intersection for the user specified cells?