[Fransje]

Dear Foamers,

I'm using the cuttingPlane class to sample information on a plane during a simulation run. I then use the .cutCells() function to retrieve the cells cut by my plane. But in the labelList returned, every cell cut is sampled twice instead of once.. Any idea why??
My code is:

point pnt_1(1.5, 1, 2);
vector direction(1, 0, 0);
plane plane_1(pnt_1, direction);

cuttingPlane cutPlane_1(plane_1, mesh);

const labelList& cutCellsPlane_1 = cutPlane_1.cutCells();

And the output of the labelList is (with Info ):

cutCellsPlane_1   face-> 1
cutCellsPlane_1   face-> 1
cutCellsPlane_1   face-> 7
cutCellsPlane_1   face-> 7
cutCellsPlane_1   face-> 13
cutCellsPlane_1   face-> 13
etc.

Is this behaviour familiar to anybody?

Kind regards,

Francois.

[olesen]

Originally Posted by Fransje

Dear Foamers,

I'm using the cuttingPlane class to sample information on a plane during a simulation run. I then use the .cutCells() function to retrieve the cells cut by my plane. But in the labelList returned, every cell cut is sampled twice instead of once.. Any idea why??
My code is:

point pnt_1(1.5, 1, 2);
vector direction(1, 0, 0);
plane plane_1(pnt_1, direction);
 
cuttingPlane cutPlane_1(plane_1, mesh);
 
const labelList& cutCellsPlane_1 = cutPlane_1.cutCells();

And the output of the labelList is (with Info ):

cutCellsPlane_1   face-> 1
cutCellsPlane_1   face-> 1
cutCellsPlane_1   face-> 7
cutCellsPlane_1   face-> 7
cutCellsPlane_1   face-> 13
cutCellsPlane_1   face-> 13
etc.

Is this behaviour familiar to anybody?

Kind regards,

Francois.

I'd guess that cuttingPlane is indeed triangulating as stated in the documentation:
http://foam.sourceforge.net/doc/Doxy....html#_details

[Fransje]

Dear Mark,

Thank you for you reply!

I did indeed read in the doxygen documentation that there could be triangulation involved, but I don't understand why that should return a .cutCells() list returning every cell twice. Especially on a Cartesian orthogonal grid..

Could you explain to me why triangulation would give me such a result?

Kind regards,

Francois.

[olesen]

Originally Posted by Fransje

Dear Mark,

Thank you for you reply!

I did indeed read in the doxygen documentation that there could be triangulation involved, but I don't understand why that should return a .cutCells() list returning every cell twice. Especially on a Cartesian orthogonal grid..

Could you explain to me why triangulation would give me such a result?

Kind regards,

Francois.

In the general case (ie, a polyhedral mesh), the cutting plane would return various weird polygon shapes. For sampling purposes (eg, calculating sums or averages) triangulated or non-triangulated doesn't matter much. For post-processing, however, it is usually much easier to handle triangles instead of weird polygon shapes.

The following code snippet from cuttingPlane.C explains what is happening inside:

 
            // the cut faces are usually quite ugly, so always triangulate
            label nTri = f.triangles(cutPoints, dynCutFaces);
            while (nTri--)
            {
                dynCutCells.append(cellI);
            }

If you don't mind hacking your own version (just to get going), you can swap the above code with the following:

const bool triangulate = true; // temporary hack for Francois
 
            // the cut faces are usually quite ugly, so optionally triangulate
            if (triangulate)
            {
                label nTri = f.triangles(cutPoints, dynCutFaces);
                while (nTri--)
                {
                    dynCutCells.append(cellI);
                }
            }
            else
            {
                dynCutFaces.append(f);
                dynCutCells.append(cellI);
            }

This means, however, that everywhere that uses cuttingPlane now doesn't triangulate. You'll have to decide yourself if this is okay for you.

The good news is that triangulation/non-triangulation of cuttingPlane is an option in the next OpenFOAM release, so you won't need this hack in later releases.

[Fransje]

Dear Mark,

Once again, thank you for your swift reply!

As I was using an orthogonal grid for the time being and needed an fix fast to get my simulation running, I wrote a simple routine filtering out the cells listed twice in my cell list. Not very elegant, I know, but it did the job until I found a better way of doing it. Now that the urgency to get my simulation running is gone, I can search for a better solution.

So thank your piece of code! It definitely looks like a more elegant solution! I will implement it, see how it works, and if I need triangulation, I will know which boolean I have to change to get it back.

Kind regards,

Francois.

[olesen]

Originally Posted by Fransje

Dear Mark,

Once again, thank you for your swift reply!

As I was using an orthogonal grid for the time being and needed an fix fast to get my simulation running, I wrote a simple routine filtering out the cells listed twice in my cell list. Not very elegant, I know, but it did the job until

For these types of quick filtering, you can just use a HashSet or a PackedBoolList to create a unique list ... shouldn't need very many lines of code.

[Fransje]

Dear Mark,

Well, it appears I keep getting more reasons to thank you!

I had only vaguely heard of HashSet until today, although it appears to be used in quite fundamental handling of objects in OpenFOAM...
But after greping and doxygen-ing around and what not, I came up indeed with a very short piece of code to do my filtering!

labelHashSet myHashList(myCutCellsLabelList);
labelList myFilteredList = myHashList.toc();

Was that the type of code you were having in mind, or is there an even cleaner way of doing the job?

Thank you for the comments!

Kind regards,

François.

[olesen]

Originally Posted by Fransje

Dear Mark,

Well, it appears I keep getting more reasons to thank you!

I had only vaguely heard of HashSet until today, although it appears to be used in quite fundamental handling of objects in OpenFOAM...
But after greping and doxygen-ing around and what not, I came up indeed with a very short piece of code to do my filtering!

labelHashSet myHashList(myCutCellsLabelList);
labelList myFilteredList = myHashList.toc();

Was that the type of code you were having in mind, or is there an even cleaner way of doing the job?

Thank you for the comments!

Kind regards,

François.

Yes Hashes are quite useful in general, not just in OpenFOAM.
You indeed have the type of code I was thinking about. For some cases the PackedBoolList might be faster (you can try the applications/test/PackedList2/ if you'd like).

The only way to shorten your code (which is already quite short) would be to skip intermediates:

const labelList myFilteredList(labelHashSet(cutPln.cutCells()).toc());

See if it works - sometimes the compiler gets a bit cranky and may not like it.

[Fransje]

Dear Mark,

Well, thank you, once more! It works!

It was a bit odd because the compiler started by being cranky, not wanting to compile the code, so I changed it back to the original code, compiled, and retried with exactly the same inline code once more, and it worked.. Dont ask..

On a totally different topic, you wouldn't happen to be familiar with the reduce() function to reconstruct data across processors by any chance?

Kind regards,

Francois.

[olesen]

Originally Posted by Fransje

...
On a totally different topic, you wouldn't happen to be familiar with the reduce() function to reconstruct data across processors by any chance?

You shouldn't have any problems finding examples from within OpenFOAM itself ("git grep --cached" is very helpful).

For your sample plane application, it could get used something like this:

   // count all faces
        label nFaces = pln.faces().size();
        reduce(nFaces, sumOp<label>());
 
   // find max value (and known a priori to be > 0)
         scalar maxU(0);
         if (pln.faces().size())
         {
some operation to get a max value from the plane
          }
          reduce(maxU, maxOp<scalar>());

Take a look in the code for many, many more examples.