CFD Online Discussion Forums - Ray Tracing in OpenFOAM

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Ray Tracing in OpenFOAM (https://www.cfd-online.com/Forums/openfoam-programming-development/127449-ray-tracing-openfoam.html)

MMC15

December 11, 2013 11:29

Ray Tracing in OpenFOAM

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	December 11, 2013 13:05

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

December 12, 2013 06:53

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	December 12, 2013 07:11

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

Code:



#include "polyMesh.H"

#include "meshTools.H"

#include "treeDataFace.H"

Code:

    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

    );

Code:

        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

December 12, 2013 08:54

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

Code:

DynamicList<pointIndexHit> pHitList;

Code:

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	December 12, 2013 09:03

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:

Code:

const point pStart = ...; 

const point pEnd = ...; 

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

MMC15

December 12, 2013 09:27

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

Code:

 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	December 12, 2013 09:36

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

Code:

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

December 12, 2013 09:42

OK! Thank you a lot for your useful replies!

Bye!

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

Yahoo

May 16, 2014 12:38

including header files

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

May 16, 2014 13:04

Hi,

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

Code:

-I$(LIB_SRC)/meshTools/lnInclude

in your options file.

Yahoo

May 27, 2014 17:04

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

alexeym

May 28, 2014 03:01

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	May 28, 2014 03:32

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:

Code:

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

February 22, 2015 06:00

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

I'm doing ray tracing, using octree:

Code:

indexedOctree<treeDataFace> faceTree     

(         

     treeDataFace(false, mesh),         

     allBb, // overall search domain         

     8, // maxLevel         

     10, // leafsize         

     3.0 // duplicity     

);

and then

Code:

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:

Code:

treeDataFace(false, mesh),

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

Code:

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

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

David*

December 8, 2015 08:21

Quote:

Originally Posted by mirx (Post 494455)

Code:

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

Code:

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

April 21, 2017 10:58

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

December 14, 2017 13:14

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

All times are GMT -4. The time now is 03:56.