[Edy]

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 :

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 :

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]

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:

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]

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]

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

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

Hi again,

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

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

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]

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]

Originally Posted by Edy

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

writeCellCentres

does what I need here.

[Sherlock_1812]

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]

Hi Srivaths,


You can extract the coordinates in one direction with

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]

Originally Posted by Sherlock_1812

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]

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

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

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

[Lieven]

Dear Yeru,

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

const vectorField& c = patch().Cf();
const scalarField xcoord(c & Vector::X);

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

forAll(coord,faceI)
{
   scalar x = xcoord[faceI];
}

afterwards

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

Cheers,

L

[Yeru]

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

 // 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!

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

does the trick.

[Lieven]

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]

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]

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

Originally Posted by Edy

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 :

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 :

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]

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]

#include "Vector.H"

[TurbJet]

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?