CFD Online Discussion Forums - How to assign value to a single cell?

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

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

- - How to assign value to a single cell? (https://www.cfd-online.com/Forums/openfoam-programming-development/182317-how-assign-value-single-cell.html)

How to assign value to a single cell?

Hi!
I'm developing custom solver based on rhoCentralFoam. There I have three main fields initialized as follows:

Code:

volScalarField& p = thermo.p();

const volScalarField& T = thermo.T(); 

volVectorField U

(

    IOobject

    (

        "U",

        runTime.timeName(),

        mesh,

        IOobject::MUST_READ,

        IOobject::AUTO_WRITE

    ),

    mesh

);

Later, I create some scalars:

Code:

      

scalar tmpT, tmpP, tmpUx, tmpUy, tmpUz;

and read them from some file.
Now I need to assign them them to cell subset, something like this:

Code:

forAll (selectedCellSet, celli)

{

  T[celli] = tmpT;  

  p[celli] = tmpP;  

  U[celli] = Vector(tmpUx, tmpUy, tmpUz);

}

This code generates various error messages, including 'cannot convert' statements and 'assignment of read-only location' (for T). What would be the correct way to assign them?

Code:

const volScalarField& T = thermo.T();

==>

Code:

volScalarField& T = thermo.T();

Thanks, Zeppo! That solved the problem with 'assignment of read-only location'.

The main question still stands: how to assign scalar value to T[celli]? Maybe somehow through dimensionedScalar / dimensionedVector construction?

Code:

volScalarField& T = thermo.T();

T.internalField() = tmpT; // for OF 3.0

//T.primitiveFieldRef() = tmpT; // for OF 4.0

Quote:

Originally Posted by Zeppo (Post 632333)

Code:

T.internalField() = tmpT; // for OF 3.0

Wouldn't it set the entire field to tmpT instead of only T[celli]?
And what is the correct syntax for U[celli], which is vector?

In addition to previous post:

Quote:

Originally Posted by Zeppo (Post 632329)

Code:

const volScalarField& T = thermo.T();

==>

Code:

volScalarField& T = thermo.T();

It turned out that this doesn't work:

Code:

createFields.H:11:29: error: binding ‘const volScalarField 

{aka const Foam::GeometricField<double, Foam::fvPatchField,

 Foam::volMesh>}’ to  reference of type ‘Foam::volScalarField& {aka

 Foam::GeometricField<double, Foam::fvPatchField, 

Foam::volMesh>&}’ discards qualifiers

Hey,

did you try it with a copy-process ?

Code:

volScalarField T = thermo.T();

T.internalField()[celli] = value;

Depending on the structure of thermo.T() you should be able to assign a scalar or a dimensionedScalar. That´s the clean solution.
Otherwise you can try to assign a non dimensioned scalar with:

Code:

T.internalField()[celli].value() = value;

That´s the dirty one :rolleyes:

You will need a volVectorField to assign a vector or a dimensionedVector to it. In addition you can assign single vector components by the functions vector.x() or vector.y() or vector.z().

Please note that this can vary depending on the OF versions you use.

Thanks, Traction!
The first way compiles OK in foam-extend 3.2 (which is based on OF 3, I believe). It generates some weird patchy field value assignments, which is probably have no relation to this thread (some errors in cell indexes, probably).

For now I'm abandoning this approach in favor of re-using some part of setFields code, which is much simpler. I've posted some problems I have with it in another thread: https://www.cfd-online.com/Forums/op...tml#post632372

Quote:

Originally Posted by OlegSutyrin (Post 632364)

In addition to previous post:

Code:

volScalarField& T = thermo.T();

It turned out that this doesn't work:

Code:

createFields.H:11:29: error: binding ‘const volScalarField 

{aka const Foam::GeometricField<double, Foam::fvPatchField,

 Foam::volMesh>}’ to  reference of type ‘Foam::volScalarField& {aka

 Foam::GeometricField<double, Foam::fvPatchField, 

Foam::volMesh>&}’ discards qualifiers

The reason for the code being not compiled might be that the thermo was declared as a const reference like

Code:

autoPtr<rhoThermo> thermoPtr = new rhoThermo(vfMesh);

const rhoThermo& thermo = thermoPtr();

If so, then

Code:

volScalarField& T = thermo.T();

won't work.

Quote:

Originally Posted by OlegSutyrin (Post 632360)

Code:

T.internalField() = tmpT; // for OF 3.0

Wouldn't it set the entire field to tmpT instead of only T[celli]?
And what is the correct syntax for U[celli], which is vector?

Code:

T.internalField() = tmpT ;

sets every cell's value to be tmpT

Quote:

Originally Posted by OlegSutyrin (Post 632360)

And what is the correct syntax for U[celli], which is vector?

Code:

U[celli] = vector(Ux, Uy, Uz);

Quote:

Originally Posted by Traction (Post 632366)

Hey,
did you try it with a copy-process ?

Code:

volScalarField T = thermo.T();

T.internalField()[celli] = value;

Code:

volScalarField T = thermo.T();

is copying from thermo.T() to T. And now you have 2 self-contained volScalarFields which are independent on each other. Modifications made to T will not propagate to thermo. This might not be what you want in some situations.

Zeppo, thanks again for such an elaborate answer!

The thing that puzzles me is that pThermo is created as non-const:

Code:

autoPtr<basicPsiThermo> pThermo

(

    basicPsiThermo::New(mesh) 

);

basicPsiThermo& thermo = pThermo();

but 'thermo.T()' still returns 'const volScalarField&'.

And I agree, the copy of my 'thermo.T' would be useless for me since I need to update 'rho', which is derived 'thermo'.

Quote:

Originally Posted by OlegSutyrin (Post 632377)

The thing that puzzles me is that pThermo is created as non-const:

Code:

autoPtr<basicPsiThermo> pThermo

(

    basicPsiThermo::New(mesh) 

);

basicPsiThermo& thermo = pThermo();

but 'thermo.T()' still returns 'const volScalarField&'.

You are right, T() can be called only for constant object:

Code:

const Foam::volScalarField& Foam::basicThermo::T() const

This thermo seems to be designed with the implication that you first initialise it with T then solve your equation for e or h then you can calculate T based on e or h.

Meanwhile, you can try this little hacking:

Code:

volScalarField& T = const_cast<volScalarField&>

(

    static_cast<const basicPsiThermo&>(thermo).T()

);