CFD Online Discussion Forums - Accessing dimensions (x, y, z) of cells

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

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

- - Accessing dimensions (x, y, z) of cells (https://www.cfd-online.com/Forums/openfoam-programming-development/235829-accessing-dimensions-x-y-z-cells.html)

CFDanielGER

April 30, 2021 07:37

Accessing dimensions (x, y, z) of cells

Hello everybody,

for a new hybrid turbulence model I want to implement I need the dimensions of each cell to calculate a new delta. I found a way of calculating the dimensions in each direction (x, y, z) in meter. See code attached:

Code:

  const faceList & ff = this->mesh().faces();

  const pointField & pp = this->mesh().points();

  volScalarField deltax = this->y_;

  volScalarField deltay = this->y_;

  volScalarField deltaz = this->y_;



  forAll(this->mesh().C(), celli)

  {

    const cell & cc = this->mesh().cells()[celli];

    labelList pLabels(cc.labels(ff));

    pointField pLocal(pLabels.size(), vector::zero);



    forAll(pLabels, pointi)

    {

      pLocal[pointi] = pp[pLabels[pointi]];

    }



    dimensionedScalar Convert ("Convert", dimensionSet(0, 1, 0, 0, 0, 0, 0), scalar(1));



    deltax = (Foam::max(pLocal & vector(1,0,0)) - Foam::min(pLocal & vector(1,0,0))) * Convert;

    deltay = (Foam::max(pLocal & vector(0,1,0)) - Foam::min(pLocal & vector(0,1,0))) * Convert;

    deltaz = (Foam::max(pLocal & vector(0,0,1)) - Foam::min(pLocal & vector(0,0,1))) * Convert;

  }

I think the code is working fine, at least it is compiling and a few iterations are running but I have to check further.

But my question is:

Is there a more effective way of computing the cell dimensions in x, y and z than using a forAll loop? I have the impression that this is rather time-consuming when using a large mesh.

Many thanks for your help!

Daniel

olesen

May 3, 2021 02:35

boundBox for the cell points?

BTW in your code example it appears that you are overwriting the entire fields each cell - probably not what you want.

For efficiency, it would make sense to store/access the values as a DimensionedField of vectors and dump all of it into a class as a MeshObject. This would make it demand-driven and you could handle mesh changes seamlessly.

CFDanielGER

May 3, 2021 03:12

Maybe I show the whole code so it becomes a bit more clear. I should also mention that I am a beginner with OpenFOAM coding and found this code in a other forum post and it seemed that this the thing I need: https://www.cfd-online.com/Forums/op...ize-x-y-z.html

My whole code for calculating the variable I need with the cell dimensions looks like this. For a non-moving mesh, the dimensions of each cell should only calculated at the beginning since they are not changing in time, but I am not sure yet how to solve this issue. Maybe write another file like the cubeRootVolDelta or the maxDelta. But what changes with each iteration is the vorticity, this should therefore be calculated at each time step.

Code:

template<class BasicTurbulenceModel>

tmp<volScalarField> SAEDES<BasicTurbulenceModel>::deltaOmega

(



) const

{

  volVectorField vorticity = fvc::curl(this->U());

  volScalarField magVorticity = mag(vorticity);

  volVectorField vorticityUnit = vorticity / magVorticity;



  const faceList & ff = this->mesh().faces();

  const pointField & pp = this->mesh().points();

  volScalarField deltax = this->y_;

  volScalarField deltay = this->y_;

  volScalarField deltaz = this->y_;



  forAll(this->mesh().C(), celli)

  {

    const cell & cc = this->mesh().cells()[celli];

    labelList pLabels(cc.labels(ff));

    pointField pLocal(pLabels.size(), vector::zero);



    forAll(pLabels, pointi)

    {

      pLocal[pointi] = pp[pLabels[pointi]];

    }



    dimensionedScalar Convert ("Convert", dimensionSet(0, 1, 0, 0, 0, 0, 0), scalar(1));



    deltax = (Foam::max(pLocal & vector(1,0,0)) - Foam::min(pLocal & vector(1,0,0))) * Convert;

    deltay = (Foam::max(pLocal & vector(0,1,0)) - Foam::min(pLocal & vector(0,1,0))) * Convert;

    deltaz = (Foam::max(pLocal & vector(0,0,1)) - Foam::min(pLocal & vector(0,0,1))) * Convert;



    //Info <<deltax<< endl;

    //Info <<deltay<< endl;

    //Info <<deltaz<< endl;

  }



  return sqrt( sqr(vorticityUnit.component(0)) * deltay * deltaz + sqr(vorticityUnit.component(1)) * deltaz * deltax + sqr(vorticityUnit.component(2)) * deltax * deltay );

}

Quote:

Originally Posted by olesen (Post 803002)

Where do you mean the entire fields will be overwritten each cell, sorry if i don't understand yet clearly. And may you please explain or describe further how your idea of storing and accessing the values look like?

Many thanks for your help!

Best regards,
Daniel

olesen

May 3, 2021 05:32

Quote:

Originally Posted by CFDanielGER (Post 803005)

If you examine the closing statements in your loop:

Code:

  forAll(this->mesh().C(), celli)

  {

...

// !!! HERE ! !!!



     deltax = (Foam::max(pLocal & vector(1,0,0)) - Foam::min(pLocal & vector(1,0,0))) * Convert;

    deltay = (Foam::max(pLocal & vector(0,1,0)) - Foam::min(pLocal & vector(0,1,0))) * Convert;

    deltaz = (Foam::max(pLocal & vector(0,0,1)) - Foam::min(pLocal & vector(0,0,1))) * Convert;

  }

You are assigning deltax, deltay etc, which are volume fields with a single value.
After exiting the loop, the entire field will simply have the value of the last cell, which you will not notice if the mesh is uniform.

olesen

May 3, 2021 05:48

You would need to verify, but my first idea was something like this (NOTE: quality not checked!!).

Code:

volScalarField deltax = this->y_;

volScalarField deltay = this->y_;

volScalarField deltaz = this->y_;



... starting here

const polyMesh& pMesh = this->mesh();

const label nCells = pMesh.nCells();





boundBox cellBb;

for (label celli = 0; celli < nCells; ++celli)

{

    const labelList& cPoints = pMesh.cellPoints(celli);

    cellBb.clear();

    cellBb.add(pMesh.points(), cPoints);



   deltax[celli] = (cellBb.max().x() - cellBb.min().x());

   deltay[celli] = (cellBb.max().y() - cellBb.min().y());

   deltaz[celli] = (cellBb.max().z() - cellBb.min().z());

}

CFDanielGER

May 3, 2021 07:52

Many thanks for your help, now I understand! You are completely right. I just tested your piece of code and it seems to work also it seems much faster, I just have to check if the dimensions are right. I think for this the easiest way is to create IOobjects, what do you think?

But if I now want to calculate the deltaOmega for each cell, I don't have to use a forAll loop, right? I can just calculate like this:

Code:

sqrt( sqr(vorticityUnit.component(0)) * deltay * deltaz + sqr(vorticityUnit.component(1)) * deltaz * deltax + sqr(vorticityUnit.component(2)) * deltax * deltay );

olesen

May 3, 2021 11:09

It's fine enough the way you now have it.
Since you have complete volfields after the loop, the rest should be OK.

You will Probably want to store deltaxy, deltayz, deltazx instead. Cheap to calculate within the loop, and saves extra multiplications later.

After that can revise how you create your fields. Probably want them with the class scope. For that, autoPtr is the usual container.

saicharan662000@gmail.com

May 20, 2022 07:53

Hi Daniel,

what does this->Y in your code mean??
Thanks in advance

CFDanielGER

May 20, 2022 09:45

Quote:

Originally Posted by saicharan662000@gmail.com (Post 828346)

Hi Daniel,

what does this->Y in your code mean??
Thanks in advance

That should be the wall distance if I remember correctly. It was just for initialization.

saicharan662000@gmail.com

May 20, 2022 10:05

I had written a loop similar to yours for finding deltaX and wrote return statement after the loop. But when I use deltaX in return statement it is throwing an error message that deltaX is undefined

CFDanielGER

May 23, 2022 02:08

Quote:

Originally Posted by saicharan662000@gmail.com (Post 828353)

I had written a loop similar to yours for finding deltaX and wrote return statement after the loop. But when I use deltaX in return statement it is throwing an error message that deltaX is undefined

Did you define deltaX before the loop?

saicharan662000@gmail.com

May 23, 2022 02:12

Hi Daniel,
I didn't define it before loop. Because I don't know what I should assign to deltaX as u assigned this->Y. my cose looks like this

Foam::tmp<Foam::volScalarField>
Foam::phaseChangeTwoPhaseMixtures::myPhaseChange:: kappaGradT() const
{

const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
forAll(owner, facei)
{
vector delta = C[neighbour[facei]] - C[owner[facei]];
volScalarField gradAlpha = fvc::grad(limalpha1) & mesh.C();
volScalarField alphaPrim = mag(gradAlpha);
const dimensionedScalar len("len",dimLength,mag(delta));
const dimensionedScalar small("small",dimensionSet(0, 0, 0, 0, 0, 0, 0),1e-5);

}
return
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(len))*(gradAlpha/max(alphaPrim,small))
);
}

CFDanielGER

May 23, 2022 02:24

Quote:

Originally Posted by saicharan662000@gmail.com (Post 828440)

you could try to define a vector before the loop like:

Code:



vector delta(0,0,0);

saicharan662000@gmail.com

May 23, 2022 02:32

Hi Daniel,
Thanks for your suggestion,I implemented it like this now I have no eroors or warnings,is it correct way??

Foam::tmp<Foam::volScalarField>
Foam::phaseChangeTwoPhaseMixtures::myPhaseChange:: kappaGradT() const
{

const volScalarField limalpha1(min(max((alpha1()), scalar(0)), scalar(1)));
const volVectorField& U = db().lookupObject<volVectorField>("U");
const fvMesh & mesh = U.mesh();
const unallocLabelList& neighbour = mesh.neighbour();
const unallocLabelList& owner = mesh.owner();
const volVectorField& C = mesh.C();
const volScalarField& T = alpha1().db().lookupObject<volScalarField>("T");
vector delta(0,0,0);
const dimensionedScalar len("len",dimLength,0);
forAll(owner, facei)
{
vector delta = C[neighbour[facei]] - C[owner[facei]];
volScalarField gradAlpha = fvc::grad(limalpha1) & mesh.C();
volScalarField alphaPrim = mag(gradAlpha);
const dimensionedScalar len("len",dimLength,mag(delta));
const dimensionedScalar small("small",dimensionSet(0, 0, 0, 0, 0, 0, 0),1e-5);

}
return
(
(limalpha1*lambda1_+ (1-limalpha1)*lambda2_) * ((T-Tsat_)/(len))*(gradAlpha/max(alphaPrim,small))
);
}

is delta in loop getting used? Why I am doubting this is, I have set both delta and len to 0 before the loop ,will the get updated

CFDanielGER

May 23, 2022 02:54

Quote:

Originally Posted by saicharan662000@gmail.com (Post 828443)

Yes, that should be working. You are just defining delta and len beforehand so that it will be recognized for the loop. During the loop it will then get updated. You could print len or delta during your simulation so you can see if the results are meaningful. I am not sure how to print anymore but there are some posts related to this topic, just google for it.

saicharan662000@gmail.com

May 23, 2022 03:33

Hi daniel,
Thanks for helping me, And I am able to print them and I am getting correct results

CFDanielGER

May 23, 2022 05:11

Quote:

Originally Posted by saicharan662000@gmail.com (Post 828451)

Hi daniel,
Thanks for helping me, And I am able to print them and I am getting correct results

Perfect, no worries!

saicharan662000@gmail.com

May 23, 2022 05:44

Hi Daniel,
Can you look this post and suggest me what's wrong in it.
https://www.cfd-online.com/Forums/op...ing-error.html
Thanks in advance

Wang Shang

November 21, 2023 10:26

Hi Olesen,

If I want to write it as a post-processing program, is that possible?

Regards.

All times are GMT -4. The time now is 20:37.