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

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

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

Reply
 
LinkBack Thread Tools Display Modes
Old   October 15, 2010, 07:31
Default Problem / Misunderstanding with mesh.C() : get the cell centers
  #1
Edy
Member
 
Join Date: Sep 2010
Posts: 35
Rep Power: 6
Edy is on a distinguished road
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
Edy is offline   Reply With Quote

Old   October 16, 2010, 21:03
Default
  #2
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   October 18, 2010, 08:25
Default
  #3
Edy
Member
 
Join Date: Sep 2010
Posts: 35
Rep Power: 6
Edy is on a distinguished road
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
Edy is offline   Reply With Quote

Old   October 18, 2010, 08:43
Default
  #4
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   October 18, 2010, 08:53
Default
  #5
Edy
Member
 
Join Date: Sep 2010
Posts: 35
Rep Power: 6
Edy is on a distinguished road
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());
    }
Edy is offline   Reply With Quote

Old   October 18, 2010, 09:07
Default
  #6
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 14
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   October 18, 2010, 09:48
Default
  #7
Edy
Member
 
Join Date: Sep 2010
Posts: 35
Rep Power: 6
Edy is on a distinguished road
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
Edy is offline   Reply With Quote

Old   July 9, 2012, 08:05
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 Edy View Post
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.

Last edited by iamed18; July 9, 2012 at 08:45. Reason: Resolved
iamed18 is offline   Reply With Quote

Old   August 1, 2013, 12:51
Default writing to internalField
  #9
Senior Member
 
Srivathsan N
Join Date: Jan 2013
Location: India
Posts: 101
Rep Power: 4
Sherlock_1812 is on a distinguished road
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?
__________________
Regards,

Srivaths
Sherlock_1812 is offline   Reply With Quote

Old   August 2, 2013, 02:10
Default
  #10
Senior Member
 
Lieven
Join Date: Dec 2011
Location: Mol, Belgium
Posts: 295
Rep Power: 13
Lieven will become famous soon enough
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
Lieven is offline   Reply With Quote

Old   August 2, 2013, 04:35
Default
  #11
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
Quote:
Originally Posted by Sherlock_1812 View Post
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?
Search for funkySetFields
Bernhard is offline   Reply With Quote

Old   June 30, 2014, 10:10
Default Alternating scalar value on boundaryField
  #12
New Member
 
Dominik Pöltl
Join Date: Jul 2013
Location: Hamburg
Posts: 14
Rep Power: 4
Yeru is on a distinguished road
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?
Yeru is offline   Reply With Quote

Old   June 30, 2014, 10:19
Default
  #13
Senior Member
 
Lieven
Join Date: Dec 2011
Location: Mol, Belgium
Posts: 295
Rep Power: 13
Lieven will become famous soon enough
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
Lieven is offline   Reply With Quote

Old   June 30, 2014, 10:37
Default
  #14
New Member
 
Dominik Pöltl
Join Date: Jul 2013
Location: Hamburg
Posts: 14
Rep Power: 4
Yeru is on a distinguished road
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.

Last edited by Yeru; July 3, 2014 at 09:43. Reason: Solution found and added
Yeru is offline   Reply With Quote

Old   June 30, 2014, 13:52
Default
  #15
Senior Member
 
Lieven
Join Date: Dec 2011
Location: Mol, Belgium
Posts: 295
Rep Power: 13
Lieven will become famous soon enough
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
Lieven is offline   Reply With Quote

Old   July 1, 2014, 04:50
Default Saving coordinates in a seperate file
  #16
New Member
 
Dominik Pöltl
Join Date: Jul 2013
Location: Hamburg
Posts: 14
Rep Power: 4
Yeru is on a distinguished road
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
Yeru is offline   Reply With Quote

Old   April 3, 2015, 12:23
Post
  #17
Member
 
Muhammad Usman
Join Date: Feb 2014
Posts: 67
Rep Power: 3
13msmemusman is on a distinguished road
Hye can you describe how will you define Ua????? i cant understand it
Quote:
Originally Posted by Edy View Post
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
13msmemusman is offline   Reply With Quote

Old   May 7, 2015, 16:42
Default
  #18
New Member
 
rafath
Join Date: Jun 2014
Location: mumbai
Posts: 17
Rep Power: 3
R_21 is on a distinguished road
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.
R_21 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
Cell Reynolds Number laliong Main CFD Forum 11 March 19, 2015 12:42
Adiabatic and Rotating wall (Convection problem) ParodDav CFX 5 April 29, 2007 19:13
Cell face values computation un unstructured grids Sergio Rossi Main CFD Forum 2 May 28, 2006 10:04
May I use global coordinate of cell centers in subroutine Posdat H. Shen CD-adapco 4 March 29, 2000 09:20
extremely simple problem... can you solve it properly? Mikhail Main CFD Forum 40 September 9, 1999 09:11


All times are GMT -4. The time now is 15:01.