[gwierink]

Dear all,

I am working with a project where we would like to implement a surface capturing method that requires information of the coordinate system in the form x(i,j,k). I have understood that OpenFOAM uses a cell list of the form x(i), so you do not really know the ijk position of a cell. The grid we want to use is a simple orthogonal grid. I would like to know the following:

1. Does anyone know of a method to retrieve the "ijk information" of each cell in OpenFOAM?
2. How can the cell number (or cell position info) be accessed in OpenFOAM? Is it a "real cell number" or a e.g. a pointer to a location in memory?

Many thanks in advance!

[ngj]

Hi

Well, you can achieve the information from the following building blocks:

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

forAll( mesh.C(), celli)
{
    // 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)
    {
         // Using own and nei you can determine  x(i,j,k)
    }

}

What method is it you are considering to implement?

Best regards,

Niels

[gwierink]

Hi Niels,

Fantastic! Thank you for your quick reply!

The method I want to implement is called the Constrained Interpolation Profile (CIP) method. It is basically a subgrid reconstruction method which you can use to capture an interface. I am interested in modelling interaction between gas-liquid interfaces and particles (or drops of a third phase/species). By using CIP on an orthogonal grid I hope to be able to implement more coupled physics easier, like Marangoni effects, Eh-pH, electrophoretic forces etc.

Apologies for the noob-question , but would you have some hints on how to determine x(i,j,k) from own and nei? Thanks in advance!

[ngj]

Hi Gijs

Okay, that is interesting. I did an implementation of a DDR scheme following Harvie & Fletcher earlier this year. It has some better properties with respect to keeper a clean interface, however my implementation unfortunately seems to be somewhat faulty.

With respect to neighbouring/owner, then

// Cell labels
label owner = own[cProp[facei]];
label neigh = nei[cProp[facei]];

if (owner == celli)
    // then the neighbour is actaully neigh
else // neigh == celli
  // then the neighbour is actually owner

Hope it makes sense.

Best regards,

Niels

[gwierink]

Hi Niels,

Thanks for your answer. Due to some other issues (related to the physico-chemistry we want to solve) we decided to use only part of the CIP method, so we don't really need the x(i,j,k) coordinates at this point. Sorry for the hassle about that .

To implement the (partial) CIP method I want to try transporting a solid using interFoam first, as a test. I intend to kind of freeze the solid and liwuid around it and then translate it according to the flow field. This is to prevent diffussion of the solid (which will otherwise happen with this VOF method). In order to do this I need to browse through all the faces and check whether there is on one cell center water and on the other side air. If I then find the faces that are "interface faces" I want to "grow" a (few) cell layer(s) around the solid object and then transport the "grown solid block". Then, I want to integrate the forces on the faces of the grown interface. Would you have some hints how to implement this? I have some pseudo-code below, but lots is still missing ...

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


// Loop through all cell faces
forAll(mesh.Sf(), facei)
{ 
    // Check for all faces whether it has on one side
    // a cell center with alpha1 = 0 and a cell center
    // with alpha1 = 1 on the other side.
    forAll( mesh.C(), celli)
    {
        // Grow a cell layer in the direction of alpha1 = 0
    }  
}

Thanks in advance!

[ngj]

Hi Gijs

Unfortunately I cannot help you with your problem, however if the object you are going to consider remain a simple geometry throughout the computation, then you might want to consider moving the boundary between the two fluids, e.g. [1].

Z. Tukovic has further developed a FAM-routine, where you solve given differential equations on a surface. As you have chemistry in your problem it might be relevant so forth some of the reactions are limited to take place on the interface.

Good luck,

Niels

[1]: Tukovic, Z. and Jasak, H. (2007). Updated Lagrangian finite volume solver
for large deformation dynamic response of elastic body. Transactions of
FAMENA, 31(1), 55–70.

[gwierink]

Hi Niels,

Thanks for your help.

Z. Tukovic has further developed a FAM-routine, where you solve given differential equations on a surface. As you have chemistry in your problem it might be relevant so forth some of the reactions are limited to take place on the interface.

Indeed, good idea. I will keep digging .

[ngj]

You are welcome. I remembered having a link to the mentioned article, see [1]. Especially section 5 might interest you.

Further, regarding FAM-framework it is already implemented in 1.5-dev and there is a tutorial where the dispersion of a species on a curved surface is modelled.

Best regards,

Niels


[1] http://powerlab.fsb.hr/ped/kturbo/Op...tionFAMENA.pdf

[gwierink]

Hi Niels,

Is there an easy way to loop through all (internal) faces in a mesh? Something like the following perhaps (although this doesn't compile)?

const unallocLabelList& listOfFaces = mesh.faces();
forAll(listOfFaces, faceI)
{
    // do stuff
}

[ngj]

Hi Gijs

I have looped over the internal faces of a surfaceScalarField, e.g. phi, when doing this. I guess you need to flux anyway for your two phase solver, so that might be a way to go.

Best regards,

Niels

[gwierink]

Hi Niels,

Right, so you're suggesting to interpolate the volScalarField onto the faces and then loop through them? I mean, phi is already the face flux, isn't it? I am dealing with cell centered values (alpha). I guess it can be done, but is it efficient?

[ngj]

What I mean is that when I have needed what you are doing, I was considering the direction of the flux over the cell faces, hence I needed to consider the sign of phi on every face relative to the orientation of the hence. This made it a obvious choice.

Could you elaborate on what it is you need?

Best regards,

Niels

[gwierink]

Sure, sorry for the fuzzyness .

I have a phase field alpha1, which is set to 1 inside a small blob and 0 in the main fluid. What I want to do is "pad" the blob with another phase, which is actually the same as the main fluid, just name differently for convenience. Let's call that phase beta. So, I want to loop through all faces and compare alpha1 on each side. When alpha1 is significantly different on one side than the other, it is an interface and I want to "pad" the adjacent main fluid as beta.

Now I have this:

forAll(alpha1, cellI)
{
    beta[cellI] = alpha1[cellI];

    label ownerCell     = mesh.faceOwner()[cellI];
    label neighbourCell = mesh.faceNeighbour()[cellI];

    // then, compare owner and neighbour
}

beta is also a volScalarField and simply a working copy of alpha1. This kind of works, but for some reason only cells below the blob are padded, not on top. So my square blob of 2x2 cells becomes a strawberry . I like strawberries, but not like this . Would you have any suggestions? Thanks in advance!

[gwierink]

Hi Niels,

I got rid of my strawberry!

flux over the cell faces

As you suggested, I used the face fluxes to loop over all the internal faces! phic is already there, so I used that to find a list of faces. Then, for all those face I compared alpha on each side. Now the blob is also padded "on the top"! A big thanks for the hint!

[ngj]

You are welcome.

Have a good weekend,

Niels

[Erik]

Hi all and Gijsber,

I would like to find the minimum and maximum x, y and z coordinates in my geometry while running my solver. I can see these values from the case's blockMeshDict file but I would like to avoid manually entering these scalar values by hand in the solver code.

I find that when NOT running in parallel the following lines works fine:

//values below somewhere in the middle of the domain
scalar minX=0, minY=0, minZ=0, maxX=0, maxY=0, maxZ=0;

forAll( mesh.cellCentres(), cellPrev)
{
const cell& cProp = mesh.cells()[cellPrev];
forAll(cProp, facei)
{
if (mesh.faceCentres()[cProp[facei]].x() < minX)
minX=mesh.faceCentres()[cProp[facei]].x();
if (mesh.faceCentres()[cProp[facei]].y() < minY)
minY=mesh.faceCentres()[cProp[facei]].y();
if (mesh.faceCentres()[cProp[facei]].z() < minZ)
minZ=mesh.faceCentres()[cProp[facei]].z();

if (mesh.faceCentres()[cProp[facei]].x() > maxX)
maxX=mesh.faceCentres()[cProp[facei]].x();
if (mesh.faceCentres()[cProp[facei]].y() > maxY)
maxY=mesh.faceCentres()[cProp[facei]].y();
if (mesh.faceCentres()[cProp[facei]].z() > maxZ)
maxZ=mesh.faceCentres()[cProp[facei]].z();

}
}

However, as I run these lines in a parallel simulation they do not yield the min and max coordinate values for the respective processor domain. I take it that this is because some processors do not own their boundary face, but a cell of the neighboring processor cell does, and therefor these bordering vertices are not seen.

Any suggestion for how to retrieve the min and max coordinate values is appreciated!

Cheers,
Erik Arlemark

[l_r_mcglashan]

If the mesh doesn't change during your simulation then the easiest thing to do is to have a separate 'preprocessing' application that works things like this out for you in serial. Otherwise you have to do something more complicated.

[Erik]

Thank you Laurence R. McGlashan,

I chose to get the coordinates by creating an extra dictionary in my case, where I put in all the user inputs being specific for my case.

Cheers,
Erik

[gascortado]

Hi all.

I'm working in developing a fluid structure interaction solver for solid-fluid. I was looking into obtaining the coordinates of each of the cells at the interface so this post helps.

While searching the forums I notice how people seem to know all kind of classes in openfoam and what they actually do. My big question is: Is there a list or guide or something where I can learn about the different classes that I can use in OF? At this point I would like to know about class that operate on the mesh.

I will appreciate any help. Thanks

[mirko]

The source code documentation (see here) has a list of all the classes.

Mirko