CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   loop through structured mesh by i,j,k (

akidess July 20, 2010 11:27

loop through structured mesh by i,j,k
Hello all,

I know I can loop through all cells of a mesh using forAll(mesh.cells(),celli). However, this won't tell me anything about the geometric location of the neighbors. I'd like to only use the neighbors on a certain axis for my further calculations. Does anyone have an idea?

In pseudo code, in 2D:

for (i =0; i < mesh.max_i; i++)
    for (j =0; j < mesh.max_j; i++)
          // calculation involving mesh.cell[j-1], mesh.cell[j], mesh.cell[j+1]

deepsterblue July 20, 2010 12:41

I suppose you could look at cell-faces, take the dot-product of the unit face-normal with the axis-direction that you want, and check its magnitude. You can then choose the owner/neighbour of those faces to get the adjacent cell.

akidess July 21, 2010 05:01

Thanks, that got me started! Now I'm facing a weird problem though: Every now and then, instead of a sensible neighbour cellID, I get something way off, eg.


inner circle, celli: 932
neigh: 962 //makes sense
inner circle, celli: 932
neigh: 1751339886 //can't be, grid only has like 1200 cells, using it to access a field naturally leads to a seg fault

Any Idea why that would be happening and how to avoid it?


vector axis(0, 1, 0);

    forAll( mesh.C(), celli)
            // Neighbour and owner cells to given face
            const unallocLabelList & nei = mesh.neighbour();
            const unallocLabelList & own = mesh.owner();

            // This is a list of faces of which the cell is build
            const cell& cProp = mesh.cells()[celli];

            // Loop over all faces
            forAll(cProp, facei)
                label owner = own[cProp[facei]];
                label neigh = nei[cProp[facei]];

                //find faces in y-direction
                if (mag(mesh.Sf()[facei] & axis) > 0)
                    Info << "celli: " << celli << endl;
                    if (owner == celli) {
                        Info << "neigh: " << neigh << endl; 
                    } else { // neigh == celli
                        Info << "owner: " << owner << endl;

ngj July 21, 2010 05:10

That would occur, if you are requesting a neighbour for a boundary face. You should make a conditional like


if (cProp[facei] < own.size())
    // Can ask for both owner and neighbour
    // Only owner exists



akidess July 21, 2010 05:27

So my guess is that the large number are boundary cells and could be just filtered out by comparing with mesh.size() [edit: Sorry Niels, didn't see your post while writing this one]. However, I now noticed that I'm not seeing two neighbors as expected on a 2D mesh, but three! Whats going on? Here is an image of cellIDs I got out of paraview.

Looking at cell 31, I want to get 1 and 61 as output. Instead I'm getting:
celli: 31
neigh: 61
owner: 31

celli: 31
owner: 31
neigh: NULL

celli: 31
neigh: 31
owner: 30

Why is there an inner face, and why am I getting the left neighbor instead of the lower one?

All times are GMT -4. The time now is 06:52.