[iamed18]

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

$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:

	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]

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:

const volVectorField& Cf = mesh.Cf();

and work with this.

The cell centres you get in the same way:

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:

    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]

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]

Originally Posted by stevenvanharen

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]

Agreed.

I think if you change this line:

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

into:

scalar distance = mag(Cf[facei] - C[owner[facei]]);

Than you have what you need?

[iamed18]

Originally Posted by tomislav_maric

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:

const volVectorField& Cf = mesh.Cf();

and work with this.

The cell centres you get in the same way:

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:

    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]

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]

Originally Posted by tomislav_maric

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]

Originally Posted by iamed18

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]

Originally Posted by tomislav_maric

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]

Originally Posted by iamed18

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]

Originally Posted by tomislav_maric

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

[tomislav_maric]

Originally Posted by iamed18

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

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

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:

   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]

Originally Posted by tomislav_maric

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 . So, when I see the following for the center-to-face distance in a few cases scattered throughout the output:

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]

@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]

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

#include "fvCFD.h"

I would assume that this header is required for some additional surface field functionality.

[saicharan662000@gmail.com]

Hi maric,
I added a subroutine that calculates cell center to center distance in x direction simliar to yours. But I am not getting the correct outPut. When I run a test case on 1D it is running fine. But for 2D its not working properly.

Can you tell what's wrong in my code.
I already posted it , I am attaching the link of my post here :capturing cell centers


Thanks in advance