CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Calculate Cell Centre to Cell Face Distances (http://www.cfd-online.com/Forums/openfoam-programming-development/104422-calculate-cell-centre-cell-face-distances.html)

iamed18 July 9, 2012 10:46

Calculate Cell Centre to Cell Face Distances
 
I'm in need of finding the distance from every cell center to the center of each of its faces. I've been able to use the utility in
Code:

$FOAM_APP/utilities/postProcessing/miscellaneous/writeCellCentres
to write out the XYZ cell center vector, but my attempts to find the cell face vectors have been unfruitful. My current attempt was merely to mimic the former:

Code:

        surfaceVectorField cf
        (
            IOobject
            (
                "cellFaces",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh.Cf()
        );

        cf.write();

This doesn't compile, leading me to believe that I'm missing something with how I'm trying to access the Cf() method of mesh. Any ideas?

tomislav_maric July 9, 2012 11:41

Cf is a field that is stored within the mesh. It is done using lazy evaluation: the field is evaluated once, meaning that there is a dynamic memory allocation taking place, so that when you call mesh.Cf() multiple times (no changes to the mesh assumed), it will not be re-calculated.

For this reason, you need only a const ref to the field:


Code:

const volVectorField& Cf = mesh.Cf();
and work with this.

The cell centres you get in the same way:

Code:

const volVectorField& C = mesh.C();
There is no need to copy(construct) these fields, since they are kept by the mesh class. Once you get the constant references (equivalence to const volVectorField* const if you come from C ), use forAll loop on the faces, get the owner of each face, and compute the distance:

Code:

    const surfaceVectorField& Cf  = mesh.Cf();
    const volVectorField& C = mesh.C();

    surfaceScalarField faceCellDist
    (
        IOobject
        (
            "faceCellDist",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar
        (
            "0",
            dimLength,
            scalar(0)
        )
    );

    const labelList& own = mesh.owner();

    forAll (own, I)
    {
        faceCellDist[I] = Foam::mag(C[own[I]] - Cf[I]);
    }

    faceCellDist.write();

That's just to get you started.

This will only run in serial and access the distance to the cell that owns the face. If you need the distance to the neighbour as well, you can store both values in a surfaceVectorField, and keep the third component of the vector untouched. This way, you have an overhead in a value per face, but you don't get into writing GeometricFIeld<Pair, ... > . If space is a constraint (doubt it), you can have a List<Pair> with both owner-cell-face, neighbour-cell-face distance stored in a Pair. The length of the List will be equal to the number of faces to the mesh.

For parallel stuff, you need to go over boundaries and access some data accross it, and then you have the same problem that I posted regarding Pstream: you will need to communicate distances across patches "at the same time".

Tomislav

stevenvanharen July 9, 2012 12:04

I would take a look in:

src » finiteVolume » interpolation » surfaceInterpolation » surfaceInterpolation» surfaceInterpolation.C

The makeWeights function is almost what you need I think.

tomislav_maric July 9, 2012 13:35

Quote:

Originally Posted by stevenvanharen (Post 370491)
I would take a look in:

src » finiteVolume » interpolation » surfaceInterpolation » surfaceInterpolation» surfaceInterpolation.C

The makeWeights function is almost what you need I think.

Yep. But they give you a single surfaceScalarFied (no P - face - N distacne or weight pair). I don't know what the actual algorithm is about, but weights_ is the thing to use if the data is needed for linear interpolation: for equidistant orthogonal hex mesh the weight will be 0.5.

stevenvanharen July 9, 2012 14:39

Agreed.

I think if you change this line:

Code:

scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
into:

Code:

scalar distance = mag(Cf[facei] - C[owner[facei]]);
Than you have what you need?

iamed18 July 10, 2012 14:58

1 Attachment(s)
Quote:

Originally Posted by tomislav_maric (Post 370485)
Cf is a field that is stored within the mesh. It is done using lazy evaluation: the field is evaluated once, meaning that there is a dynamic memory allocation taking place, so that when you call mesh.Cf() multiple times (no changes to the mesh assumed), it will not be re-calculated.

For this reason, you need only a const ref to the field:


Code:

const volVectorField& Cf = mesh.Cf();
and work with this.

The cell centres you get in the same way:

Code:

const volVectorField& C = mesh.C();
There is no need to copy(construct) these fields, since they are kept by the mesh class. Once you get the constant references (equivalence to const volVectorField* const if you come from C ), use forAll loop on the faces, get the owner of each face, and compute the distance:

Code:

    const surfaceVectorField& Cf  = mesh.Cf();
    const volVectorField& C = mesh.C();

    surfaceScalarField faceCellDist
    (
        IOobject
        (
            "faceCellDist",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar
        (
            "0",
            dimLength,
            scalar(0)
        )
    );

    const labelList& own = mesh.owner();

    forAll (own, I)
    {
        faceCellDist[I] = Foam::mag(C[own[I]] - Cf[I]);
    }

    faceCellDist.write();

That's just to get you started.

This will only run in serial and access the distance to the cell that owns the face. If you need the distance to the neighbour as well, you can store both values in a surfaceVectorField, and keep the third component of the vector untouched. This way, you have an overhead in a value per face, but you don't get into writing GeometricFIeld<Pair, ... > . If space is a constraint (doubt it), you can have a List<Pair> with both owner-cell-face, neighbour-cell-face distance stored in a Pair. The length of the List will be equal to the number of faces to the mesh.

For parallel stuff, you need to go over boundaries and access some data accross it, and then you have the same problem that I posted regarding Pstream: you will need to communicate distances across patches "at the same time".

Tomislav

Tomislav,

I like the idea you've put forth here, but unfortunately compilation failed on this as well. I've attached the error output as well as the code written thus far. It seems to complain about an incorrect usage of surfaceScalarField.

I appreciate all input!
Thanks,
~Ed

tomislav_maric July 10, 2012 15:32

1 Attachment(s)
Hi,

it compiles on my laptop, either I posted something wrong in the code or you copied wrong. :)

Anyhew, I'm using 2.1.x, but for this, it makes absolutely no difference. I have attached the application to this post..

Tomislav

iamed18 July 10, 2012 16:53

Quote:

Originally Posted by tomislav_maric (Post 370687)
Hi,

it compiles on my laptop, either I posted something wrong in the code or you copied wrong. :)

Anyhew, I'm using 2.1.x, but for this, it makes absolutely no difference. I have attached the application to this post..

Tomislav

Well, that appears to work quite nicely! Thank you!

As a final question about it, I do need to ask: does the order in which the distances are listed (from the cell center) follow the convention labeled in the OF manual?

tomislav_maric July 10, 2012 16:57

Quote:

Originally Posted by iamed18 (Post 370698)
Well, that appears to work quite nicely! Thank you!

As a final question about it, I do need to ask: does the order in which the distances are listed (from the cell center) follow the convention labeled in the OF manual?

Short answer: yep.

Long answer:

Well, this will get you the distance to the owner of the face. You still need to address the neighbour (you can safely put it in the same loop and separate the patchFace-owner distance in the loop over boundary faces), and for the patches that are coupled, you will need to do some communication. You may end up with two surfaceScalarFields this way, or a List<Pair<scalar,scalar>>.

Btw, what are you using this for?

iamed18 July 11, 2012 09:09

Quote:

Originally Posted by tomislav_maric (Post 370700)
Short answer: yep.

Long answer:

Well, this will get you the distance to the owner of the face. You still need to address the neighbour (you can safely put it in the same loop and separate the patchFace-owner distance in the loop over boundary faces), and for the patches that are coupled, you will need to do some communication. You may end up with two surfaceScalarFields this way, or a List<Pair<scalar,scalar>>.

Btw, what are you using this for?

I'm trying to convert file output from OpenFOAM to a rendering software that my employer maintains. To do so, I need to define, for each cell, the following:

cell center vector (x,y,z)
cell center-to-face for all 6 faces (the renderer can only handle hex, simplifying life on this end), and I must know which number represents which face

So, in sum, for each FV element, I need 9 inputs to reconstruct the geometry of the simulation.

Which brings me to: I understand (slightly) what you mean about needing to worry about the neighbour cells as well as the owners, but how will I be able to tell (in the list) which cell a face belongs to, and is there a standard fashion in which the owner/neighbour is decided?

tomislav_maric July 12, 2012 03:41

Quote:

Originally Posted by iamed18 (Post 370822)
I'm trying to convert file output from OpenFOAM to a rendering software that my employer maintains. To do so, I need to define, for each cell, the following:

cell center vector (x,y,z)
cell center-to-face for all 6 faces (the renderer can only handle hex, simplifying life on this end), and I must know which number represents which face

So, in sum, for each FV element, I need 9 inputs to reconstruct the geometry of the simulation.

Which brings me to: I understand (slightly) what you mean about needing to worry about the neighbour cells as well as the owners, but how will I be able to tell (in the list) which cell a face belongs to, and is there a standard fashion in which the owner/neighbour is decided?

In this case, you can loop over cells without a problem and you don't need to worry about owner-neighbour relationship. I thought you will need the distances for some implicit numerical calculations, in which case, you would need to loop over own-nei.

If I may just suggest, next time explain the actual task, this will save time ... :)

If I (finally) understood you right, you need a cell centre vector written down in a file and next to it, 6 distances (scalars) from 6 different faces.

In what order do you need those distances written? This is a structured mesh format, right? Do you need distance between the East West North South face, or something similar?

iamed18 July 12, 2012 11:30

Quote:

Originally Posted by tomislav_maric (Post 371008)
If I (finally) understood you right, you need a cell centre vector written down in a file and next to it, 6 distances (scalars) from 6 different faces.

In what order do you need those distances written? This is a structured mesh format, right? Do you need distance between the East West North South face, or something similar?

I apologize; I had never meant to waste any time, I just didn't know how to define my problem, perhaps. I appreciate the help none-the-less.

You've hit the nail on the head at this point, it sounds like. This is a structured mesh, and as far as order goes, if we use the cell center as the reference:
-x,-y,-z,+x,+y,+z

But if they don't come out of OF in that order, that's fine, I have to convert everything to binary anyway :D

tomislav_maric July 12, 2012 14:22

Quote:

Originally Posted by iamed18 (Post 371151)
I apologize; I had never meant to waste any time, I just didn't know how to define my problem, perhaps. I appreciate the help none-the-less.

You've hit the nail on the head at this point, it sounds like. This is a structured mesh, and as far as order goes, if we use the cell center as the reference:
-x,-y,-z,+x,+y,+z

But if they don't come out of OF in that order, that's fine, I have to convert everything to binary anyway :D

No need to apologize, just think it's faster to get help this way! :D

Okay, so basically, you have to re-create this connectivity for the whole mesh, since OF mesh has no idea about the axis alignment. So, here's a try:

Code:


  const cellList& cells = mesh.cells();
 
   
    const surfaceVectorField& Cf  = mesh.Cf();
    const volVectorField& C = mesh.C();

    vectorField basis(3);
    basis[0] = vector(1,0,0);
    basis[1] = vector(0,1,0);
    basis[2] = vector(0,0,1);

    List<List<scalar> > centreFaceDists (cells.size());

    forAll (cells, I)
    {
        const cell& c = cells[I];

        scalarField distances(c.size());

        forAll (c, J)
        {
            const vector& cf = Cf[c[J]];
            vector local  = cf - C[I];

            forAll (basis, K)
            {
                scalar dist = basis[K] & local;
                if (dist < 0)
                {
                    distances[K] = dist;
                    // The other projections will be 0.
                    break;
                }
                else if (dist > 0)
                {
                    distances[K + 3] = dist;
                    // The other projections will be 0.
                    break;
                }
            }
        }

        centreFaceDists[I] = distances;
    }

    //Info << centreFaceDists << endl;

    forAll (centreFaceDists, I)
    {
        Info << C[I].x() << " "
            << C[I].y() << " "
            << C[I].z() << " "
            << centreFaceDists[I][0] << " "
            << centreFaceDists[I][1] << " "
            << centreFaceDists[I][2] << " "
            << centreFaceDists[I][3] << " "
            << centreFaceDists[I][4] << " "
            << centreFaceDists[I][5] << endl;
    }

I did this very quickly, its ugly, but it compiles, check it please, I hope it helps. ;)

iamed18 July 12, 2012 15:07

Quote:

Originally Posted by tomislav_maric (Post 371171)
I did this very quickly, its ugly, but it compiles, check it please, I hope it helps. ;)

This is really cool! I'm still very early in my OF learning, and seeing creative uses of their tools is neat.

Now, the mesh I'm using for test is structured and regular; each cell should be a rectangular prism, and checkMesh reports my smallest volume to be about v_{min}\approx0.0045m^3. So, when I see the following for the center-to-face distance in a few cases scattered throughout the output:

Code:

132.431 -1.10069 8.64881 -0.1 -2.22045e-16 -0.10119 5.68434e-14 0.100694 0.10119
132.631 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 0.100694 0.10119
132.831 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 0.100694 0.10119
133.031 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 0.100694 0.10119
133.231 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 0.100694 0.10119
133.431 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 0.100694 0.10119
133.631 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 0.100694 0.10119
133.831 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 0.100694 0.10119
134.031 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 2.22045e-16 0.10119
134.231 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 2.22045e-16 0.10119
134.431 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 2.22045e-16 0.10119
134.631 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 2.22045e-16 0.10119
134.831 -1.10069 8.64881 -0.1 -0.100694 -0.10119 2.84217e-14 0.100694 0.10119
135.031 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 0.100694 0.10119
135.231 -1.10069 8.64881 -0.1 -0.100694 -0.10119 0.1 0.100694 0.10119

I become a little hesitant. I'm not learned enough yet to pick out why this error would happen, but perhaps you'll see it straight away. I used the code you gave me prior verbatim.

Thoughts?

tomislav_maric July 17, 2012 05:04

@iamed:

I'm not sure why this happens, I've tried it on simple box hex meshes and it seems to be working without a problem.

Anyway: you have a starting point now, time to do some debugging. :)

maninthemail October 10, 2012 06:09

I see this has been largely dealt with, but just in case anyone else is struggling...

Tomislav's code compiles in 2.0.x, using writeCellCentres as a template, if (and only if) I add

Code:

#include "fvCFD.h"
I would assume that this header is required for some additional surface field functionality.


All times are GMT -4. The time now is 05:10.