[Fauster]

Dear Foamers

I would like to know how to find the downwind cells of a cell. As mentionned here : owner and neighbor #4 :

A positive flux goes from "owner" to "neighbour".

.

At the moment I do the folowing thing :

forAll(T.internalField(),celli) // loop over cell
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const labelList& faces = mesh.cells()[celli];

forAll(faces,facei) // loop over faces of celli
{
if (own[facei] == celli)
{
///// Is nei[facei] a downwind cell of celli ?
Pout << "phi[facei] = " << phi[facei] << endl ;
// Problem : phi[facei] can be positive or negative.
}
}
}

If we assume that a "positive flux goes from "owner" to "neighbour" " I don't understant why phi[facei] can be positive or negative in a face if own[facei] == celli ?
So I suppose I am doing wrong something. Has anyone an idea how to solve this problem ?

Another question :
what is the difference between

mesh.neighbour()[facei] AND mesh.faceNeighbour()[facei]

?

Huge thanks

[frobaux]

Hi,
I think there is a problem with your indices, you should do own[faces[facei]] instead of faces[facei] if you want the owner cell of the facei-th face of the cell. Same applies to Phi[facei]=> Phi[faces[facei]]. (otherwise facei goes from 1 to not much - like 10 max)



I think that a positive flux goes from owner to neighbor just means that, when the flux is positive, the flow effectively goes from owner to neighbour. When it is negative, the flow is reversed. Basically, that just means that the normal of a face is oriented outwards (from its owner to its neighbour)


for the faceNeighbour() vs neighbour() I don't know where you find the "neighbour()" function.

faceNeighbour()[f] is returning the neighbouring cell of the face of label f (again you should use "label f=faces[facei]").



I hope I'm clear,
Fabien

[Fauster]

Dear Frobaux,

Thank you for you message. I am going to use your solution and see if it's ok. Since my first message I found out that in isoadvector there is a function wich seems to calculate downwindFaces :

void Foam::isoAdvection::setDownwindFaces
(
const label celli,
DynamicLabelList& downwindFaces
) const
{
DebugInFunction << endl;

// Get necessary mesh data and cell information
const labelList& own = mesh_.faceOwner();
const cellList& cells = mesh_.cells();
const cell& c = cells[celli];

downwindFaces.clear();

// Check all faces of the cell
forAll(c, fi)
{
// Get face and corresponding flux
const label facei = c[fi];
const scalar phi = faceValue(phi_, facei);

if (own[facei] == celli)
{
if (phi > 10*SMALL)
{
downwindFaces.append(facei);
}
}
else if (phi < -10*SMALL)
{
downwindFaces.append(facei);
}
}

downwindFaces.shrink();
}

I understand that this function fill a list with the downwind faces of a cell. If I look at the way of looping over cell and then faces I don't see so much difference with my code. What am I missing ?

For example in the code of linearUpwind :

https://github.com/OpenFOAM/OpenFOAM...linearUpwind.C

The build of the surfaceScalarField returned by the "correction" function is obtained by looping over all faces of the mesh. From my understanding there is no [faces[facei]]. Clearly I am missig something. Could you go more in details please ?

for the faceNeighbour() vs neighbour() I don't know where you find the "neighbour()" function.

It seems that neighbour() is for example here :

https://github.com/OpenFOAM/OpenFOAM...linearUpwind.C

Huge thank for you answer.
P

[frobaux]

Yes, you are right faceNeighbour() is defined in polyMesh.H and neighbour() is defined in fvMesh.H. I don't know the difference, but as fvMesh is a child of polyMesh, you have access to both of them..
and the same is true for owner() and faceOwner() by the way. However be careful, the type they return are different, faceOwner returns a list of face when owner returns a list of labels.

The difference between your code and the downwind cell is the indexing:
you do

forAll(c, fi) // fi goes from 0->8 if 8 faces on cell c


        const scalar phi = faceValue(phi_, fi); // whatever the cell, you wil get only phi[0] .. phi[8] so the 8 first flux of the faces of mesh 

            if (own[fi] == celli)
    {
         ....

    }

}

when they convert to global index: (line facei=c[fi])

forAll(c, fi)

         const label facei = c[fi];   // for the, lets say 3rd face of the cell c, c[3] is the global index of that face! 
        const scalar phi = faceValue(phi_, facei);

            if (own[facei] == celli)
    {
         ....

    }

}