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/)
-   -   Problem / Misunderstanding with mesh.C() : get the cell centers (http://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)
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

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)
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


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.


All times are GMT -4. The time now is 08:50.