# Distance of a 'North' and 'South' surface of a cell for 2D mesh

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

October 6, 2020, 10:22
Distance of a 'North' and 'South' surface of a cell for 2D mesh
#1
Senior Member

René Thibault
Join Date: Dec 2019
Posts: 112
Rep Power: 5
Hi everyone!

My problem could be easy for some people, but I personally don't have any idea on how to code this thing here.

I would like to calculate, for a 2D mesh, the distance 'z' from the wall to the 'north' (Zn) and 'South' (Zs) faces of a cell (see image in attachment).

Im able to obtain the distance for a center of a cell (ZN, ZP, ZS) with 'mesh_.C' like so:
Code:
`const vectorField& CellC = mesh_.C();`
But I don't have any clue on how to get the distance 'z' for the 'North' and 'South' faces of the same cell. What I have understand so far, and I'm maybe wrong, is that I need to find from a particular cell the indices what are the 'North' and 'South' surface are compose of and calculate the normal vector in order to to calculate the height 'z' after.

But I don't have any idea whatsoever how to code this thing.

Can I have some example and some help?

Thank you everyone!
Attached Images
 Cells.png (2.7 KB, 14 views)

Last edited by Tibo99; October 6, 2020 at 14:52.

 November 27, 2020, 05:58 #2 New Member   Join Date: Feb 2019 Posts: 23 Rep Power: 6 Hi, You can loop over all cell faces and take into account only the faces which have the z component of its normal vector different from 0. Then, see which face is the N and which is S. Finally, just take the distance from the wall (assuming z=0) to the face centers. The code could be something like (for cell ID = 10, for example): Code: ```const cellList& cellFaces = mesh.cells(); const surfaceVectorField& Sf = mesh.Sf(); const surfaceVectorField& Cf = mesh.Cf(); const surfaceVectorField& C = mesh.C(); scalar zWall = 0.0; label cellI = 10; scalar zn = 0.0; scalar zs = 0.0; forAll(cellFaces[cellI],faceI) { label faceIG = cellFaces[cellI][faceI]; scalar zSf = Sf[faceIG][2]; if (zSf != 0.0) { if (Cf[faceIG][2] > C[cellI][2]) { zn = Cf[faceIG][2] - zWall; } if (Cf[faceIG][2] < C[cellI][2]) { zs = Cf[faceIG][2] - zWall; } } }``` Notice that I haven't taken into account boundary faces. If so, the code would be slightly different.

 November 27, 2020, 09:06 #3 Senior Member   René Thibault Join Date: Dec 2019 Location: Canada Posts: 112 Rep Power: 5 Thank you very much for replying! The fact that the code don't take into account any boundary faces in the domain is not that a big deal for me since I use it, for the moment, for a 2D channel. My 1st question is, because I need to apply this from the 2nd cell to the 4th cells from the wall in the entire domain, how I could modify the loop to take this into account? The 2nd one is, is this option work with parallel computing? I was thinking using 'topoSet' to get the ID for these cells and then, I would use your loop inside a new loop in order to apply it to all of the cells that 'topoSet' has selected. If this approach is right, I don't know at this point how to get all of these cells ID, stored them in a constant ''IDList'' (for instance) and using it after in the new loop. I keep working on it. Even though, my main issue is pretty much solved with the loop you suggested. Thank you very much again and be safe. Regards, Last edited by Tibo99; November 29, 2020 at 19:32.

 November 28, 2020, 06:04 #4 Senior Member   Mark Olesen Join Date: Mar 2009 Location: https://olesenm.github.io/ Posts: 1,554 Rep Power: 36 Your use case is fairly specific, but could still be worth generalizing. Starting with a simple struct with four label members (north, south, east, west) that will hold the cellId of the respective neighbours. You could also just use a FixedList of 4 labels for this purpose. Create a List of these structs that is nCells long, with each element initialized to -1. Code: ```enum compass { NORTH = 0, SOUTH, EAST, WEST }; typedef FixedList compassType; List compassAddr(mesh.nCells(), compassType(-1));``` Next you want to loop across each mesh face. Internal faces have own/neighbour cells, boundary faces only have an owner cell. From the face you get its normal. Since you only care about the general direction, just use the area normal (possibly slightly cheaper). You don't need the scalar (dot) product, but just find the cmptMax (the max component) of the face normal. With this you know if this is east-west or north-south. The combination of the sign of that component and the owner/neighbour relationship tells you which sub elements of the respective compassAddr entries to set. Now that the addressing is complete, you can now handle everything else that you need. A simple example (untested) Code: ```// cellCentres() for polyMesh const pointField& cc = mesh.C(); scalar dist = 0; label celli = someCellID; // walk four to the east for (label i=0; i < 4; ++i) { const label next = compassAddr[celli][compass::EAST]: if (next == -1) break; dist += mag(cc[next] - cc[celli]); celli = next; }``` I think this is enough for you to work with.

 November 28, 2020, 06:16 #5 Senior Member   Mark Olesen Join Date: Mar 2009 Location: https://olesenm.github.io/ Posts: 1,554 Rep Power: 36 Just re-read your first post and noticed that you actually wanted the centre to face distances. In this case you need to store the face id instead of the cell id in the north/south/east/west slots. This will give you an immediate calculation of centre to face, but when walking to the next cell you will need to use the mesh face owner/neighbour information for the connectivity and to walk in the correct direction. This is probably more flexible than what I originally posted.

 November 29, 2020, 13:29 #6 Senior Member   René Thibault Join Date: Dec 2019 Location: Canada Posts: 112 Rep Power: 5 Hi Olesen, I maybe should have mentioned in the first post that I still considered myself at a beginner level with OF. Regarding the fact that the solution you suggested is generalized and probably the best one for my application, I have to admit, I don't exactly understand the code. I feel that this is an another level for me, which I'm not there yet. Sorry about that. But, here is a thread that I posted and it's related to. That should help everyone to understand why I asked how to code this. Applying correction to the k-e transport equation Thank you very much for your time and be safe! Last edited by Tibo99; November 30, 2020 at 10:57.