CFD Online Discussion Forums - Problem / Misunderstanding with mesh.C() : get the cell centers

Page 1 of 2

Show 40 post(s) from this thread on one page

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Problem / Misunderstanding with mesh.C() : get the cell centers (https://www.cfd-online.com/Forums/openfoam-programming-development/81084-problem-misunderstanding-mesh-c-get-cell-centers.html)

Edy	October 15, 2010 07:31

Problem / Misunderstanding with mesh.C() : get the cell centers

Hej,

This is probably very simple but I experienced unexpected results with the .C() function.

It is supposed to return the cell centers of the mesh it is applied to. I think that I am right so far...

However i have a 2D geometry defined for x between 0.0 and 0.0127 and y between 0.0 and 1.7, as you can see from my blockMeshDict below :

Code:

vertices

(

    (0      0   0)//0

    (0.0127 0   0)//1

    (0.0127 1.7 0)//2

    (0      1.7 0)//3

    (0      0   0.1)//4

    (0.0127 0   0.1)//5

    (0.0127 1.7 0.1)//6

    (0      1.7 0.1)//7



);



blocks

(

    hex (0 1 2 3 4 5 6 7) (20 170 1) simpleGrading (1 1 1)

);

Therefore if I consider a cell center, its x-value should be between 0.0+delta and 0.0127+delta, simply because the cell has a certain dimension and the cell center cannot be located at the boundary.

However, i wrote this piece of code :

Code:

volVectorField centres = Ua.mesh().C();

volScalarField x = centres.component(0);

          Info<< min(x) << nl << endl;

          Info<< max(x) << nl << endl;

In my terminal i get:
min(x)=0.0
max(x)=0.0127

So basically it means that there are cell centers located at the boundaries, which is nonsense!

Could you please help me and tell me what i have not properly understood, and how to solve it.

Thanks in advance.
Best

/Ed

santiagomarquezd

October 16, 2010 21:03

Hi, mesh.C() gives you a slicedVolVectorField which has internalField and boundaryField, the values in the internalField are the cell centres and in the boundaryField you have the face centres of boundary patches, so the min value of this field will be the x coordinate of a face centre. The code in fvMeshGeometry.C reads:

Code:

00115 void fvMesh::makeC() const

00116 {

00117     if (debug)

00118     {

00119         Info<< "void fvMesh::makeC() : "

00120             << "assembling cell centres"

00121             << endl;

00122     }

00123 

00124     // It is an error to attempt to recalculate

00125     // if the pointer is already set

00126     if (CPtr_)

00127     {

00128         FatalErrorIn("fvMesh::makeC()")

00129             << "cell centres already exist"

00130             << abort(FatalError);

00131     }

00132 

00133     CPtr_ = new slicedVolVectorField

00134     (

00135         IOobject

00136         (

00137             "C",

00138             pointsInstance(),

00139             meshSubDir,

00140             *this,

00141             IOobject::NO_READ,

00142             IOobject::NO_WRITE,

00143             false

00144         ),

00145         *this,

00146         dimLength,

00147         cellCentres(),

00148         faceCentres()

00149     );

00150 

00151 

00152     // Need to correct for cyclics transformation since absolute quantity.

00153     // Ok on processor patches since hold opposite cell centre (no

00154     // transformation)

00155     slicedVolVectorField& C = *CPtr_;

00156 

00157     forAll(C.boundaryField(), patchi)

00158     {

00159         if (isA<cyclicFvPatchVectorField>(C.boundaryField()[patchi]))

00160         {

00161             // Note: cyclic is not slice but proper field

00162             C.boundaryField()[patchi] == static_cast<const vectorField&>

00163             (

00164                 static_cast<const List<vector>&>

00165                 (

00166                     boundary_[patchi].patchSlice(faceCentres())

00167                 )

00168             );

00169         }

00170     }

00171 }

note in lines 157-170 when boundaryField is set.

Regards.

Edy	October 18, 2010 08:25

Hi Santiago,

Thanks for your answer. I understand what the problem is.
But do you have any trick to solve this? Because I would like to obtain a field of the cell centers without the values at the boundaries, and then use this field and multiply it by other parameters to obtain a volVectorField that is used in my UEqns.

Thanks in advance!

/Edouard

santiagomarquezd

October 18, 2010 08:43

Hi Edouard, I was doing a similar work some months ago and used:

Code:

// Cell centroid coordinates

const volVectorField& centres = Ua.mesh.C().internalField();

since mesh.C() gives a you a kind of volVectorField you can extract its internalField in order to obtain only cell centres in volume mesh.

Regards.

Edy	October 18, 2010 08:53

Hi again,

Thanks again Santiago. I have solved the problem changing the x-values at the boundary, as can be seen below.

Code:

    // Take the cell centres

    volVectorField centres = Ua.mesh().C();



    volScalarField x = centres.component(0);



    // Correction of the x-field to avoid x=0 at the wall boundary

    forAll(x.boundaryField(),patchi)

    {

        x.boundaryField()[patchi] = 0.5*(max(x).value()+min(x).value());

    }

santiagomarquezd

October 18, 2010 09:07

In this case you are losing info, I think is more clear to extract the internalField because FOAM already gives you what are looking for.

Regards.

Edy	October 18, 2010 09:48

Hi,

You're right. I lose info even if I dont need it. Your solution is "cleaner", i will use it.

Thanks again for your very quick and useful answers!
Best,

/Edouard

iamed18

July 9, 2012 08:05

Quote:

Originally Posted by Edy (Post 279633)

Hi,

You're right. I lose info even if I dont need it. Your solution is "cleaner", i will use it.

Thanks again for your very quick and useful answers!
Best,

/Edouard

I know it's several years later, but this is the most relevant place I could think to ask: how can I get this mesh.C().internalField() to be saved to file?

EDIT: Alas, the shell command

Code:

writeCellCentres

does what I need here.

Sherlock_1812

August 1, 2013 12:51

writing to internalField

Hi all,

Just like reading values of the internal field, is it possible to write (assign) linearly increasing values of a scalar to the internalField on all the cell centres along one coordinate direction?

Lieven

August 2, 2013 02:10

Hi Srivaths,

You can extract the coordinates in one direction with

Code:

scalarField Xcoord = mesh.C().component(Vector::X);

You can do the same for the Y and Z coordinates. Now this field can be used to set your scalar values.

Cheers,

L

Bernhard

August 2, 2013 04:35

Quote:

Originally Posted by Sherlock_1812 (Post 443349)

Search for funkySetFields

Yeru	June 30, 2014 10:10

Alternating scalar value on boundaryField

Dear all, especially dear Lieven,

I'm currently looking for something similar like Srivathsan.
For simulating a super focus mixer, I need to assign a scalar value to my boundary faces depending on their x and z coordinate.

Since the BC is assignes locally face by face in my 0/T-file, I can only access one face's coordinates at a time.
So the question is simply:

When the expression I use is evaluated face by face, can I access the face center coordinates via

Quote:

Xcoord = mesh.C().component(Vector::X);

EDIT: No, I still get errors. Any idea?

Lieven

June 30, 2014 10:19

Dear Yeru,

You can access the x-coordinates (and similarly y and z) of a patch through:

Code:

const vectorField& c = patch().Cf();

const scalarField xcoord(c & Vector::X);

If you need the coordinates one by one, use something like

Code:

forAll(coord,faceI)

{

   scalar x = xcoord[faceI];

}

afterwards

I think this should do the trick. I haven't checked it...

Cheers,

L

Yeru	June 30, 2014 10:37

Dear Lieven,

at first, huge thanks for this express answer.
I'm afraid, it doesn't apply to my problem.

I described it in this thread as well.

The moment I want to access the coordinates, I'm "inside" the 0/T-file

Code:

 // 0/T-file

focus_inlet

{

 type           groovyBC;

variables    "circ=<value>;d_2=<value>;n=5;Xcoord=mesh.x();Zcoord=mesh.z();";

valueExpression "(<f(Xcoord,Zcoord)>    ?    1    :    0.1";

value         uniform 0.1;

So as soon as I start my simulation via scalarTransportFoam, I want each face to get assigned a value according to it's coordinates.

Your code would globally define a matrix with the face indices and their x/y-values, right?
The problem is that I can't access those via my T-file.

Confusing? I hope, I could describe the problem properly.
Dominik

EDIT: SOLVED!

Code:

variables "circ=<value>;d_2=<value>;n=5;Xcoord=pos().x;Zcoord=pos().z;";

does the trick.

Lieven

June 30, 2014 13:52

Dear Dominik,

My apologies, indeed the method I describe above can be used when adding a new BC to the sources, not by directly writing it in the e.g. 0/T-file. I missed this small detail ;-)

I have to disappoint you however, I'm not familiar with groovyBC so I won't be able to help you any further. But I'm sure there are people at this forum who can do so.

L

Yeru	July 1, 2014 04:50

Saving coordinates in a seperate file

Hello again,
I gave it some more thoughts. Lieven, I have an idea how to still use your code and then use it for the BC.

What if:
-I write a script that prints every boundary face
(index x y z)
in a n-4-matrix
-> I add a fifth column with the calculated value for the BC
-> I load this fifth column into my 0/T-file and thereby neglect implementing the calculation in groovyBC

Since I'm fairly new to OF, my questions are:
1. Where should this file be located? /constant/polyMesh or 0/T?
2. Can I easily load the fifth column via a command of by saving it as .csv?
I hope that there are some threads about the 2nd question already.

Anyway, feel invited to tell me what you think

13msmemusman

April 3, 2015 12:23

Hye can you describe how will you define Ua????? i cant understand it

Quote:

Originally Posted by Edy (Post 279344)

Code:

vertices

(

    (0      0   0)//0

    (0.0127 0   0)//1

    (0.0127 1.7 0)//2

    (0      1.7 0)//3

    (0      0   0.1)//4

    (0.0127 0   0.1)//5

    (0.0127 1.7 0.1)//6

    (0      1.7 0.1)//7



);



blocks

(

    hex (0 1 2 3 4 5 6 7) (20 170 1) simpleGrading (1 1 1)

);

Code:

volVectorField centres = Ua.mesh().C();

volScalarField x = centres.component(0);

          Info<< min(x) << nl << endl;

          Info<< max(x) << nl << endl;

R_21	May 7, 2015 16:42

Hi,

Is there any Header files to be included to use "Vector::X" because I have been getting the error,
error: ‘template<class Cmpt> class Foam::Vector’ used without template parameters.

inginheiro

June 23, 2016 06:13

#include "Vector.H"

TurbJet

January 23, 2018 14:36

how you declare mesh in other scope?

I am trying to access mesh for the sentence "forAll(mesh.C(), celli)" in other scope (another source file), but the compiler keeps telling me I have not declare mesh in the scope.

How am I gonna fix this?

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

Page 1 of 2

Show 40 post(s) from this thread on one page