# How to get cell indices and manipulate field data

 Register Blogs Members List Search Today's Posts Mark Forums Read

 June 23, 2019, 10:17 How to get cell indices and manipulate field data #1 New Member   Christian Join Date: Jun 2019 Posts: 13 Rep Power: 2 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, 05:00 #2 New Member   Christian Join Date: Jun 2019 Posts: 13 Rep Power: 2 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 ‘[Foam::label {aka int}]’ for array subscript if (T[idx] > TR) ^ qREqn.H:7:34: error: invalid types ‘[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 ‘[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 ‘[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, 08:12 #3 Super Moderator     Tobias Holzmann Join Date: Oct 2010 Location: Augsburg Posts: 2,347 Blog Entries: 6 Rep Power: 40 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, 12:36 #4 New Member   Christian Join Date: Jun 2019 Posts: 13 Rep Power: 2 Hey Tobi, thanks for your reply. I made some corrections, this code is working now: Code: ```volScalarField& T = const_cast (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.

 June 26, 2019, 04:24 #5 New Member   Christian Join Date: Jun 2019 Posts: 13 Rep Power: 2 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 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 frobaux likes this.

September 25, 2019, 11:17
#6
Senior Member

Join Date: Jan 2012
Posts: 191
Rep Power: 10
// 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 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