How to get cell indices and manipulate field data

ch1 · June 23, 2019, 11:17

Hey everybody,
the following routine written in MATLAB illustrates a basic field manipulation problem, which I am trying to implement in OF:

Code:

 1    %% Initiliaze some example fields
 2    A = rand(4,4,4); %initiliaze a scalar field with random values (mesh1)
 3    B = zeros(2,4,4); %initiliaze an empty scalar field with zeros (mesh2)
 4    x = linspace(0.25,1,4) - 0.125; %get the x-coordinates of the cell centers.
 5    %% The following needs to be implemented in OpenFOAM
 6    idx = x < 0.5; %find the indices of the cells, which are within the desired     range (x < 0.5)
 7    Aidx = A(idx,:,:); %get the corresponding field values of the cells, which are in the desired range
 8    idx2 = Aidx > 0.1; %find the indices of the cells, whose field values are higher than 0.1
 9    B(idx2) = Aidx(idx2) .* 2; %do some operation with these values and write them to the new field in mesh2
10    B(~idx2) = 0; %set the other values of that new field in mesh2 to zero

The first one is similar to the problem of mapping e.g. a volScalarField between two meshes. This is something I was also searching on the forum but I still didn't figure out how to implement it. However, to keep things simple I want to solve my problem differently: Consider two meshes A and B with different dimensions but with the same cell size (so that B could be a sub-mesh of A). I want to "copy" fields between those meshes just by appropriate addressing of the cells, see lines 6 and 7.

After that, the next step would be to do some manipulation with the field data, see lines 8 to 10.

If the MATLAB Code is unclear I can explain in more detail. I just found that it's the most compact way to describe my problem.

Best regards
Christian

ch1 · June 24, 2019, 06:00

Here is an example of how I am trying to manipulate the field:

Code:

dimensionedScalar TR ("TR", dimensionSet(0,0,0,1,0,0,0), 1000);

forAll (T,idx)
{
    if (T[idx] > TR)
    {
        qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR);
    }
    else
    {
        qR[idx] = 0.0;
    }
}

Which gives me the following compilation error:

qREqn.H: In function ‘int main(int, char**)’: /home/of/foam/foam-extend-4.0/src/foam/lnInclude/UList.H:410:36: error: overloaded function with no contextual type information for (Foam::label i=0; i<(list).size(); i++) ^ qREqn.H:3:1: note: in expansion of macro ‘forAll’ forAll (T,idx) ^~~~~~ In file included from rhoEmFoam.C:161:0: qREqn.H:5:14: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript if (T[idx] > TR) ^ qREqn.H:7:34: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); ^ qREqn.H:7:58: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); ^ qREqn.H:7:74: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR);

Tobi · June 24, 2019, 09:12

Hi,

of which kind of class does T belong? I guess it is the temperature field.

Code:

const scalarField& Tcell = T.internalField();

forAll(Tcell, celli)
{

}

By the way. You cannot compare a scalar with a dimensionedScalar which is the error you get. Either you define your limit as a normal scalar or you do extract the value of the dimensionedScalar:

Code:

    if (T[idx] > TR.value())
    {
        qR[idx] = 5600.0 * (T[idx] - TR.value()) + 181.0 * (T[idx] - TR.value()) * (T[idx] - TR.value());
    }

ch1 · June 24, 2019, 13:36

Hey Tobi,
thanks for your reply. I made some corrections, this code is working now:

Code:

volScalarField& T = const_cast<volScalarField&> (thermo.T()); //get the temperature field

//const scalar TR = 1000.0;
scalar TR(readScalar(physicalProperties.lookup("TR"))); 

const scalar c1 = 5600.0;
const scalar c2 = 181.0;

forAll(T, idx)
{
    const scalar& cellT = T[idx];
    if (cellT > TR)
    {
        qR[idx] = c1 * (cellT - TR) + c2 * sqr(cellT - TR);
    }
    else
    {
        qR[idx] = 0.0;
    }
}

The first problem was that the temperature field is by default not directly available after the enthalpy Equation is solved (I am using a modified rhoPimpleFoam solver). Now, I found how to access it, see the first line above. Also, I am converting now the cell values of the dimensionedScalar temperature field to scalar values, see the first line inside the loop, in order to avoid the illegal comparison, which you mentioned.

The last task I need to solve is the following: Find all cells of field A in meshA, which are in the range of e.g. (x < xmax), and copy the field values to a field B in meshB. Or, since I know the sizes of the meshes it might be even easier to do:
Get all values of field A but reduced by e.g. 50 cells in x-direction while keeping the number of cellValues in y and z direction. Then copy the result to the fieldB in meshB, which has the matching dimensions.

PS: I assume that you can see my second post because of your moderator status. However, I think it needs to be set to public. At least I cannot see it myself.

ch1 · June 26, 2019, 05:24

Hello,
I solved my problem by using the meshToMesh.H functions, similar to the code described in runTime mapFields from fineMesh to coarseMesh

However, I had to modify the code because I am using foam-extend 4.0.
Here is the working code:

Code:

// read the mapFieldsDict located in /case/system
IOdictionary mapFieldsDict
(
    IOobject
    (
        "mapFieldsDict",
        runTime.system(),
        runTime,
        IOobject::MUST_READ,
        IOobject::NO_WRITE,
        false
    )
);

HashTable<word> patchMap;
wordList cuttingPatches;
mapFieldsDict.lookup("patchMap") >> patchMap;
mapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh(); 
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();     

// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
    meshSource,
    meshTarget,
    patchMap,
    cuttingPatches
);

Info<< "Target mesh size: " << meshTarget.nCells() << endl;

Info<< nl
    << "Mapping fields for time " << meshSource.time().timeName()
    << nl << endl;

// specifiy source and target field 
Info<< "    interpolating qJ"      << endl;

volScalarField fieldSource(qJ);

volScalarField fieldTarget
(
    IOobject 
    (
        "qJ",
        runTime.timeName(),//meshTarget.time().timeName(),
        meshTarget,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    meshTarget,
    dimensionedScalar("qJ",dimensionSet(1,-1,-3,0,0,0,0), 0.0)
);

// Interpolate field
meshToMeshInterp.interpolate
(
    fieldTarget,
    fieldSource,
    meshToMesh::INTERPOLATE
);

// Write field
 fieldTarget.write();

Also, make sure that

Code:

#include "meshTools.H"
#include "meshToMesh.H"

is included and linked to the solver.

--------------------------------------------------

Best regards
Christian

itsme_kit · September 25, 2019, 12:17

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh();
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();

Hi Christian

I'm wondering how to define meshSource and meshTarget into emMesh and mesh respectively.

Thanks.

Quote:

Originally Posted by ch1

Hello,
I solved my problem by using the meshToMesh.H functions, similar to the code described in runTime mapFields from fineMesh to coarseMesh

However, I had to modify the code because I am using foam-extend 4.0.
Here is the working code:

Code:

// read the mapFieldsDict located in /case/system
IOdictionary mapFieldsDict
(
    IOobject
    (
        "mapFieldsDict",
        runTime.system(),
        runTime,
        IOobject::MUST_READ,
        IOobject::NO_WRITE,
        false
    )
);

HashTable<word> patchMap;
wordList cuttingPatches;
mapFieldsDict.lookup("patchMap") >> patchMap;
mapFieldsDict.lookup("cuttingPatches") >>  cuttingPatches;

// specify source and target mesh for the mapping
const fvMesh& meshSource = emMesh; //meshToMeshInterp.fromMesh(); 
const fvMesh& meshTarget = mesh; //meshToMeshInterp.toMesh();     

// Create the interpolation scheme
meshToMesh meshToMeshInterp
(
    meshSource,
    meshTarget,
    patchMap,
    cuttingPatches
);

Info<< "Target mesh size: " << meshTarget.nCells() << endl;

Info<< nl
    << "Mapping fields for time " << meshSource.time().timeName()
    << nl << endl;

// specifiy source and target field 
Info<< "    interpolating qJ"      << endl;

volScalarField fieldSource(qJ);

volScalarField fieldTarget
(
    IOobject 
    (
        "qJ",
        runTime.timeName(),//meshTarget.time().timeName(),
        meshTarget,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    meshTarget,
    dimensionedScalar("qJ",dimensionSet(1,-1,-3,0,0,0,0), 0.0)
);

// Interpolate field
meshToMeshInterp.interpolate
(
    fieldTarget,
    fieldSource,
    meshToMesh::INTERPOLATE
);

// Write field
 fieldTarget.write();

Also, make sure that

Code:

#include "meshTools.H"
#include "meshToMesh.H"

is included and linked to the solver.

--------------------------------------------------

Best regards
Christian

June 23, 2019, 11:17	How to get cell indices and manipulate field data	#1
ch1 New Member Christian Join Date: Jun 2019 Posts: 14 Rep Power: 7	Hey everybody, the following routine written in MATLAB illustrates a basic field manipulation problem, which I am trying to implement in OF: Code: 1 %% Initiliaze some example fields 2 A = rand(4,4,4); %initiliaze a scalar field with random values (mesh1) 3 B = zeros(2,4,4); %initiliaze an empty scalar field with zeros (mesh2) 4 x = linspace(0.25,1,4) - 0.125; %get the x-coordinates of the cell centers. 5 %% The following needs to be implemented in OpenFOAM 6 idx = x < 0.5; %find the indices of the cells, which are within the desired range (x < 0.5) 7 Aidx = A(idx,:,:); %get the corresponding field values of the cells, which are in the desired range 8 idx2 = Aidx > 0.1; %find the indices of the cells, whose field values are higher than 0.1 9 B(idx2) = Aidx(idx2) .* 2; %do some operation with these values and write them to the new field in mesh2 10 B(~idx2) = 0; %set the other values of that new field in mesh2 to zero The first one is similar to the problem of mapping e.g. a volScalarField between two meshes. This is something I was also searching on the forum but I still didn't figure out how to implement it. However, to keep things simple I want to solve my problem differently: Consider two meshes A and B with different dimensions but with the same cell size (so that B could be a sub-mesh of A). I want to "copy" fields between those meshes just by appropriate addressing of the cells, see lines 6 and 7. After that, the next step would be to do some manipulation with the field data, see lines 8 to 10. If the MATLAB Code is unclear I can explain in more detail. I just found that it's the most compact way to describe my problem. Best regards Christian

June 24, 2019, 06:00		#2
ch1 New Member Christian Join Date: Jun 2019 Posts: 14 Rep Power: 7	Here is an example of how I am trying to manipulate the field: Code: dimensionedScalar TR ("TR", dimensionSet(0,0,0,1,0,0,0), 1000); forAll (T,idx) { if (T[idx] > TR) { qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); } else { qR[idx] = 0.0; } } Which gives me the following compilation error: qREqn.H: In function ‘int main(int, char)’: /home/of/foam/foam-extend-4.0/src/foam/lnInclude/UList.H:410:36: error: overloaded function with no contextual type information for (Foam::label i=0; i<(list).size(); i++) ^ qREqn.H:3:1: note: in expansion of macro ‘forAll’ forAll (T,idx) ^~~~~~ In file included from rhoEmFoam.C:161:0: qREqn.H:5:14: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript if (T[idx] > TR) ^ qREqn.H:7:34: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]*’ for array subscript qR[idx] = 5600.0 (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); ^ qREqn.H:7:58: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR); ^ qREqn.H:7:74: error: invalid types ‘<unresolved overloaded function type>[Foam::label {aka int}]’ for array subscript qR[idx] = 5600.0 * (T[idx] - TR) + 181.0 * (T[idx] - TR) * (T[idx] - TR);

June 24, 2019, 09:12		#3
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Hi, of which kind of class does T belong? I guess it is the temperature field. Code: const scalarField& Tcell = T.internalField(); forAll(Tcell, celli) { } By the way. You cannot compare a scalar with a dimensionedScalar which is the error you get. Either you define your limit as a normal scalar or you do extract the value of the dimensionedScalar: Code: if (T[idx] > TR.value()) { qR[idx] = 5600.0 * (T[idx] - TR.value()) + 181.0 * (T[idx] - TR.value()) * (T[idx] - TR.value()); } __________________ Keep foaming, Tobias Holzmann

June 24, 2019, 13:36		#4
ch1 New Member Christian Join Date: Jun 2019 Posts: 14 Rep Power: 7	Hey Tobi, thanks for your reply. I made some corrections, this code is working now: Code: volScalarField& T = const_cast<volScalarField&> (thermo.T()); //get the temperature field //const scalar TR = 1000.0; scalar TR(readScalar(physicalProperties.lookup("TR"))); const scalar c1 = 5600.0; const scalar c2 = 181.0; forAll(T, idx) { const scalar& cellT = T[idx]; if (cellT > TR) { qR[idx] = c1 * (cellT - TR) + c2 * sqr(cellT - TR); } else { qR[idx] = 0.0; } } The first problem was that the temperature field is by default not directly available after the enthalpy Equation is solved (I am using a modified rhoPimpleFoam solver). Now, I found how to access it, see the first line above. Also, I am converting now the cell values of the dimensionedScalar temperature field to scalar values, see the first line inside the loop, in order to avoid the illegal comparison, which you mentioned. The last task I need to solve is the following: Find all cells of field A in meshA, which are in the range of e.g. (x < xmax), and copy the field values to a field B in meshB. Or, since I know the sizes of the meshes it might be even easier to do: Get all values of field A but reduced by e.g. 50 cells in x-direction while keeping the number of cellValues in y and z direction. Then copy the result to the fieldB in meshB, which has the matching dimensions. PS: I assume that you can see my second post because of your moderator status. However, I think it needs to be set to public. At least I cannot see it myself.