[Friederike]

Hi FOAMers,
I need to loop over all faces of a cell and collect fluxes. I have a 2D case so i set the boundaries with face normal in y-direction to empty.

I use the following code:

PHP Code:

// loop over all cells:
forAll(mesh.C(), cellI) 
{
   Info << "******* CellID: " << cellI << "*******"<< endl;

  //Getting list of all faces of current cell
  const labelList& faces = mesh.cells()[cellI];

  //loop over all faces of current cell
  forAll( faces, faceI )
  {   
    if (mesh.isInternalFace(faces[faceI]))
     {   
       Info << "internal faceI: " << faceI << "    mesh.Sf()[faceI]: " << -mesh.Sf()[faceI] << "    mesh.magSf()[faceI]: " << mesh.magSf()[faceI] << endl;
     }   
     else
     {   
       Info << "boundary faceI: " << faceI << "    mesh.Sf()[faceI]: " << -mesh.Sf()[faceI] << "    mesh.magSf()[faceI]: " << mesh.magSf()[faceI] << endl;
     }   

  } //move on to next face  
  Info << " " << endl;
}//move on to next cell 


The output of this is:

PHP Code:

******* CellID: 1*******
internal faceI: 0    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
internal faceI: 1    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1
boundary faceI: 2    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
boundary faceI: 3    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1
boundary faceI: 4    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
internal faceI: 5    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1

******* CellID: 2*******
internal faceI: 0    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
internal faceI: 1    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1
boundary faceI: 2    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
boundary faceI: 3    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1
boundary faceI: 4    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
internal faceI: 5    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1

******* CellID: 3*******
internal faceI: 0    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
internal faceI: 1    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1
boundary faceI: 2    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
boundary faceI: 3    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1
boundary faceI: 4    mesh.Sf()[faceI]: (-1 -0 -0)    mesh.magSf()[faceI]: 1
internal faceI: 5    mesh.Sf()[faceI]: (-0 -0 -1)    mesh.magSf()[faceI]: 1 


My question now are:
Why are all vector components negative (it stays like that for the whole grid)?
Why are there no vectors with a y-Component? Is this because of the pseudo 2D grid?
I get the same output with a zeroGradient boundary condition for the third dimension.
I don't want to sum over the boundaries which are set to empty. I guess I could also use this as a condition but I still would like to know why mesh.Sf() behaves so wired.

I have looked for some time now but could not find an answer...

Thanks for any help and suggestions!!!
Cheers
friederike

[stater]

Hi Friederike,
I know it is an old thread, but actually i am stuck on the same issues, so I am interested to know if you have already found the answers to your questions
Thank you

[Jerryfan]

Hi Friederike,


Are you still stuck with this problem? I found at least two problems with the above code:
First of all, The index of

mesh.Sf()[faceI]

shouldn't be faceI. It should be faces[faceI] so that you can get the global face index.
Secondly, you cannot get the surface area vector for boundary faces as Sf() only contains area vectors for internal faces.

As for why all the vector components stay negative, it is because of the mesh (all the surface normal directions are positive). For example, face 0 seems to be negative for cell 1 in x direction, but because of the owner and neighbor rule, it's always pointing from cell 0 (owner) to cell 1 (neighbor) such that it's still (-1, 0, 0).

As for why y-component is always 0, the reason is that you will not be able to get the surface vector based on Sf() even if you correct the index as suggested above. You have to access the boundary faces through the fvBoundaryMesh variable boundary_.

I hope it can help you.

[stater]

Hi,

For the first question. I think the negative value of vector components is due to the presence of the "-" character in the code just near the mesh.Sf()[faceI]

Best regards

[Jerryfan]

Hi stater,


You are right. As I mentioned above, all the surface normal directions are positive. After putting am "-" sign in front, all becomes negative. In addition to this, what I was trying to say above is that some surfaces "seem" to have negative direction, but because of the owner-neighbour convention, they are still positive.

[ashish.svm]

How to get area vector for boundary faces.

As mentioned in the first question in this thread, If a particular face is an internal face then only mesh.Sf() works.

[Zeppo]

Originally Posted by ashish.svm

How to get area vector for boundary faces.

const fvBoundaryMesh& boundary = mesh.boundary();
forAll(boundary, i)
{
    const fvPatch& patch = boundary[i];
    const vectorField& patchSf = patch.Sf();
}

or

const volVectorField::GeometricBoundaryField& boundaryField = mesh().Sf().boundaryField();
forAll(boundaryField, i)
{
    const fvPatchVectorField& patchField = boundaryField[i];
    const vectorField& patchSf = patchField.internalField();
}

[ashish.svm]

Hi Sergie,

Thanks for your reply. I tried the following code

forAll( cellIfaces, i ) // loop over all faces in the cellIndex
{
if(mesh.isInternalFace(cellIfaces[i])) // If face is internalFace
{
vector faceNormal = mesh.Sf()[cellIfaces[i]] / mesh.magSf()[cellIfaces[i]];
}
else // If face is boundaryFace
{
label patchID = mesh.boundaryMesh().whichPatch(cellIfaces[i]);
label faceID = mesh.boundaryMesh()[patchID].whichFace(cellIfaces[i]);
vector faceNormal = mesh.Sf().boundaryField()[patchID][faceID] / mesh.magSf().boundaryField()[patchID][faceID];
}
}


This compiles successfully, but during case run, it shows an error because one of the patches was defined as empty patch type. I am getting the following error:

--> FOAM FATAL ERROR:
attempt to access element from zero sized list

Can you please suggest some way resolve or get around with this type of error.