CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Calculate Cell Centre to Cell Face Distances

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   July 9, 2012, 10:46
Default Calculate Cell Centre to Cell Face Distances
  #1
Member
 
Edward Leonard
Join Date: May 2012
Location: Calumet, MI
Posts: 31
Rep Power: 5
iamed18 is on a distinguished road
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?
iamed18 is offline   Reply With Quote

Old   July 9, 2012, 11:41
Default
  #2
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
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_maric is offline   Reply With Quote

Old   July 9, 2012, 12:04
Default
  #3
Senior Member
 
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 6
stevenvanharen is on a distinguished road
I would take a look in:

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

The makeWeights function is almost what you need I think.
stevenvanharen is offline   Reply With Quote

Old   July 9, 2012, 13:35
Default
  #4
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
Quote:
Originally Posted by stevenvanharen View Post
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.
tomislav_maric is offline   Reply With Quote

Old   July 9, 2012, 14:39
Default
  #5
Senior Member
 
Steven van Haren
Join Date: Aug 2010
Location: The Netherlands
Posts: 149
Rep Power: 6
stevenvanharen is on a distinguished road
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?
stevenvanharen is offline   Reply With Quote

Old   July 10, 2012, 14:58
Default
  #6
Member
 
Edward Leonard
Join Date: May 2012
Location: Calumet, MI
Posts: 31
Rep Power: 5
iamed18 is on a distinguished road
Quote:
Originally Posted by tomislav_maric View Post
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
Attached Files
File Type: gz cellData.tar.gz (2.9 KB, 17 views)
iamed18 is offline   Reply With Quote

Old   July 10, 2012, 15:32
Default
  #7
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
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
Attached Files
File Type: zip myCellCentreFaceCentreDist.zip (47.9 KB, 41 views)
tomislav_maric is offline   Reply With Quote

Old   July 10, 2012, 16:53
Default
  #8
Member
 
Edward Leonard
Join Date: May 2012
Location: Calumet, MI
Posts: 31
Rep Power: 5
iamed18 is on a distinguished road
Quote:
Originally Posted by tomislav_maric View Post
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?
iamed18 is offline   Reply With Quote

Old   July 10, 2012, 16:57
Default
  #9
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
Quote:
Originally Posted by iamed18 View Post
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?
tomislav_maric is offline   Reply With Quote

Old   July 11, 2012, 09:09
Default
  #10
Member
 
Edward Leonard
Join Date: May 2012
Location: Calumet, MI
Posts: 31
Rep Power: 5
iamed18 is on a distinguished road
Quote:
Originally Posted by tomislav_maric View Post
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?
iamed18 is offline   Reply With Quote

Old   July 12, 2012, 03:41
Default
  #11
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
Quote:
Originally Posted by iamed18 View Post
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?
tomislav_maric is offline   Reply With Quote

Old   July 12, 2012, 11:30
Default
  #12
Member
 
Edward Leonard
Join Date: May 2012
Location: Calumet, MI
Posts: 31
Rep Power: 5
iamed18 is on a distinguished road
Quote:
Originally Posted by tomislav_maric View Post
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
iamed18 is offline   Reply With Quote

Old   July 12, 2012, 14:22
Default
  #13
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
Quote:
Originally Posted by iamed18 View Post
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:

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.
tomislav_maric is offline   Reply With Quote

Old   July 12, 2012, 15:07
Default
  #14
Member
 
Edward Leonard
Join Date: May 2012
Location: Calumet, MI
Posts: 31
Rep Power: 5
iamed18 is on a distinguished road
Quote:
Originally Posted by tomislav_maric View Post
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?
iamed18 is offline   Reply With Quote

Old   July 17, 2012, 05:04
Default
  #15
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 259
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
@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.
tomislav_maric is offline   Reply With Quote

Old   October 10, 2012, 06:09
Default
  #16
New Member
 
Bill
Join Date: Jun 2011
Location: UK
Posts: 14
Rep Power: 6
maninthemail is on a distinguished road
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.
maninthemail is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Access cell face center / cell vertices lichmaster OpenFOAM Programming & Development 7 May 31, 2014 02:31
Hydrostatic Pressure and Gravity miliante OpenFOAM Running, Solving & CFD 132 October 7, 2012 22:50
How to calculate the cell area Le FLUENT 0 February 18, 2007 23:15
Cell face values computation un unstructured grids Sergio Rossi Main CFD Forum 2 May 28, 2006 10:04
calculate cell volume, center...? Paul Main CFD Forum 5 June 21, 2003 12:55


All times are GMT -4. The time now is 19:27.