[jackdoubleyou]

Hi,


I'm creating a new solver, based on pisoFoam, that requires me to perform spatial averaging on a series of x-z planes. I've successfully created the required planes, and created cuttingPlanes from all of them.



The mesh is created using blockMesh, and has 180 points in the x direction and 40 in the z direction, so each plane should intersect 7200 cells. When I use cuttingPlane.cutCells(), instead of returning the cells that have been cut by the plane, the first 7200 cells from the mesh are given. The relevant code is given below.

yLocs is a 1D std::vector containing all of the y values I want to average, and I've confirmed that these are correct. The plane normals and origins are also being correctly used. "cells" is a labelList as I don't want to average across my whole domain, only a part of it.

labelList cutPlanes[yLocs.size()];

int i = 0;

for (const auto &yValue: yLocs)
{
  point origin(-0.018,yValue,0.002);

  vector nVec(0,1,0);

  plane nplane(origin, nVec);

  cuttingPlane plane1(nplane, mesh, 0, cells);

  Info << plane1.refPoint() << endl;

  const labelList& cellsOnPlane = plane1.cutCells();

  /*forAll(cellsOnPlane, celli)
  {
    Info << celli << endl;
  }*/

  cutPlanes[i] = cellsOnPlane;

  i++;

 }

I've been through this thread (Write cells and data intersecting a plane cuttingPlane) and not been able to fix my problem. I'm using OpenFOAM v5.x.


Thanks in advance,
Jack

[raumpolizei]

Hey,
have you tried foamToVTK + viewing the cellids in paraview? You could then check at least what should be the output. That may help reason on your implementation.
Also, I'm not sure but shouldn't cutPlanes be some kind of container<labelList> ? I don't understand what this line is supposed to do:

labelList cutPlanes[yLocs.size()];

Seems like you are mixing access operator[] and initialization which may lead to undefined behaviour.

std::vector<labelList> cutPlanes(yLocs.size());

Is this what you originally wanted to express?
Best regards
RP

[jackdoubleyou]

Hi,


I haven't tried that, will give it a go.


The std::vector line is what I originally tried, however it would not compile. I think it was because there is no implementation of std::vector that knows how to deal with labelLists.


labelList cutPlanes[yLocs.size()]


is supposed to initialise an empty array, of length yLocs.size(), of labelLists. It appears to have worked, although I agree it is not a great solution.


Thanks,
Jack

[raumpolizei]

That's strange, for me this works fine (v16.12)

#include "fvCFD.H"
#include<vector>
int main(int argc, char *argv[])
{
    std::vector<labelList> test(4);

    for( labelList& ilist: test)
    {
        ilist = labelList{1,2,3};
    }

    return 0;
}

Ah yeah, so much into C++ right now that I almost forgot about raw arrays . This should also be fine.

[jackdoubleyou]

Yeah you're right, I can't reproduce the error I was previously having. It looks like either method works, but will go with std::vector, seems neater.


Still gives me the same problem though. I've loaded in the vtk, and the cell IDs that are being returned are not the correct ones. Someone in the lab has suggested that .cutCells() returns the cell IDs local to the plane, which will always be 0-7199, instead of the cell ID in the global mesh. I can't find anything online suggesting this should be the case however.


Jack

[jackdoubleyou]

I've fixed the problem, it was to do with how I was referencing the cell list.


Instead of

forAll(cells, celli)
{
    code
}

it needed an additional line

forAll(cells, celli)
{
    const label cell = cells[celli]

    code
}

I was using the index of the labelList, instead of the cell ID contained within it.

[raumpolizei]

Riight, this forAll macro is so...unintuitive . I think this one is the only reason I, from time to time, have to lookup ranged based for loops. Glad you solved it.

[olesen]

Originally Posted by jackdoubleyou

Hi,
I'm creating a new solver, based on pisoFoam, that requires me to perform spatial averaging on a series of x-z planes. I've successfully created the required planes, and created cuttingPlanes from all of them.

If you are going to do this cutting operation a lot (and it sounds like you are), you will probably want to reconsider what you are doing.

The cuttingPlane finds intersection cells *and* attempts to create a face loop for these intersections. From your description, it doesn't really sound like you need faces. You just want to have a list of intersected cells. In fact, you don't actually really want to have a list of the cells, instead you just want to know which cells are intersected to define that plane and be able to loop over them.

There is some similar code put into OpenFOAM-v1812 that you may want to take a look at and do something similar. I looked at this recently, and the cutting problem is much easier. (Reliably constructing faces from the intersections, handling degenerate cells etc, isn't so reliable). The other reason to avoid cuttingPlane is that older versions (ie, older than 1812) worked by calculating mesh edge cuts, which means mesh edges are created for the entire mesh - this can be an expensive operation.

I think that all you want to know is which cells are intersected by the plane. This is what is happening in src/sampling/surface/cutting/cuttingPlaneCuts.C


The comments there explain the stages

// Classify sides of plane (0=BACK, 1=ONPLANE, 2=FRONT) for each point

This is saved as 2-bits per mesh point (PackedList).

Loop over each mesh face (once):

// Check for face/plane intersection based on crossings
// Took (-1,0,+1) from plane::sign and packed as (0,1,2).
// Now use for left shift to obtain (1,2,4).
//
// Test accumulated value for an intersection with the plane.

If the face intersects the plane, you mark that cell as being 'cut'. If a cell has more than one face cut, the plane intersects the cell. You could be strict and insistent on three or more cuts, but there is no benefit and just more overhead.

For count the number of times a cell cut has happened, don't waste memory actually counting the number of cuts. Instead use two separate bitSets to count your state. The reason for the bitSet is that when this little algorithm is done we have a bitSet as its output.

I really like bitSet. For these types of things it is really good for memory (1-bit per location) and gives you fast checks for testing if any bits are set, population count to get how many etc. For your case, you'll want to use the iterator directly.
Eg

bitSet cellsCut(some_algorithm_like_above);

scalar avgVal = 0;


for (const label cellId : cellsCur){
      avgVal += someField[cellId];
}

const label nTotalCut = returnReduce(cellsCuts.count(), sumOp<label>());

if (nTotalCut)
{
     reduce(sumOp<scalar>());
     avgVal /= nTotalCut;
}

This should give you good performance, low overhead and provide the result that you need. Another benefit of this approach (see the cuttingPlane calcCellCuts method in 1812), is that you can also define an input bitSet to apply additional restrictions on the intersection. Eg, if you only want plane/mesh intersections that are also within some other geometrical region (eg, bounding box) or mesh region (eg zones).


If it makes sense, could also think about breaking this private method out.


I hope this gives you some new ideas.


Cheers,
/mark

[olesen]

Originally Posted by jackdoubleyou

it needed an additional line

forAll(cells, celli)
{
    const label cell = cells[celli]
}

I was using the index of the labelList, instead of the cell ID contained within it.

For standard containers, where you don't need the local indexing number you should prefer direct iteration. This avoids a method call 'size()' during every iteration and you don't need to use an indirection to get at your values. I don't say that this will speed up your code, but makes it easier to read (I think). Thus your code,

for (const label celli : cells)
{
    Info<< "using cell: " << celli << nl;
}

You can use exactly the same trick for HashSet too. But you will need to be careful trying this on a HashTable since those iterators are already what you normally need and the dereferenced '*iter' returns the value, whereas the stl versions return a key/value pair. In most of these cases, you should then be using the `forAllConstIters()` macro, which has two parameters.


/mark