[akidess]

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]

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]

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]

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
}
else
{
    // Only owner exists
}

Cheers,

Niels

[akidess]

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?