
[Sponsors] 
How to get cell indices and manipulate field data 

LinkBack  Thread Tools  Search this Thread  Display Modes 
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 xcoordinates 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 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; } } qREqn.H: In function ‘int main(int, char**)’: /home/of/foam/foamextend4.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, 08:12 

#3 
Super Moderator

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) { } 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<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 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 xdirection 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 foamextend 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(); Code:
#include "meshTools.H" #include "meshToMesh.H"  Best regards Christian 

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:


Thread Tools  Search this Thread 
Display Modes  

