Find the neighbour cell in specific direction

[Kevin]

Hi there!
Is there a way to find the neighbor of a cell in a desired direction in cs? Example: i have cell #1 and i wanna find it's neighbour cell in z or y direction.
I found cs_geom_closest_point() but has no options for selecting in the desired direction.
Thanks

[Yvan Fournier]

Hello,

For this you need to loop over the neighboring cells, access their coordinates, and compare them with those of the current cells.

This is less practical than for a structured code, but in a purely unstructured code, there is no implict alignment of the mesh on a given direction.

In C, you can use members of cs_glob_mesh_adjacencies (see Doxygen documentation) to loop directly on neighboring cells. In Fortran you only have the face->cells connectivity so it is less practical.

Best regards,

Yvan

[Kevin]

Thank you mr.Fournier,
I just want to calculate DX+ & DZ+which needs the distances between two neighbour cells.
Is this available in cs or i shoud use cs_glob_mesh_adjacencies() and some loops?
For example, how does the cs calculate CFL number itself which needs the DX?

[Yvan Fournier]

Hello,

I am not sure how we compute the CFL (it is in the dtttvar.f90 routine), but it seems we build a Matrix based on the finite volume discretization for this, so the distances appear in FV operators, not directly as a finite difference.

Distances between neighboring cells are used in several places (such as in gradient computations).

Best regards,

Yvan

[Kevin]

Thanks admin.
I have a mesh which is generated by extruding a 2d grid in z direction with 20 layers with equal spacing in z direction. So it has cubic cells and also is homogeneous in that direction.


I found cs_glob_mesh->i_face_cells which gives the two cells that are adjacent by their common interior faces.
So, as my cells are cubic, each cell has 6 interior faces. So, if I loop over all of the interior faces I can find the neighbor cells in specific direction.
As the mesh is homogeneous in z direction, if z coordinates of two adjacent cells are equal, then they are in the same xy plane. So to test this i wrote a code to set the value of a field named Z equal to the z coordinate of the cell. So the adjacent cells on the same xy plane should have equal value of Z.

Code: Select all

 const cs_lnum_2_t *i_face_cells = (const cs_lnum_2_t *)cs_glob_mesh->i_face_cells;
        cs_lnum_t n_faces;
        cs_lnum_t cell_id1, cell_id2, face_id;
        const cs_lnum_t n_i_faces = cs_glob_mesh->n_i_faces;
        const cs_real_3_t *restrict cell_cen = (const cs_real_3_t *restrict)cs_glob_mesh_quantities->cell_cen;
        cs_real_t *Z = (cs_real_t *)cs_field_by_name("Z")->val;

  for (face_id = 0; face_id < n_i_faces; face_id++) {
    cell_id1 = i_face_cells[face_id][0];                                                                          
    cell_id2 = i_face_cells[face_id][1];                                                                          
    if((cell_cen[cell_id1][2]==cell_cen[cell_id2][2]))
      {
      Z[cell_id1] = cell_cen[cell_id1][2];
      Z[cell_id2] = cell_cen[cell_id2][2];
      }
      }

The code looks to be ok but gives zero value for some cells as you can see in the image it seems that loop is skipping some cells. This is a slice generated by paraview. The blue cells have zero value while the red cells have the correct values.

Screenshot from 2020-12-05 18-25-12.png

(2.05 KiB) Not downloaded yet

I wrote another code which does the same by using loop over all cells but this one gives correct values for all cells and all of the cells are red in paraview without any blue cells.

Code: Select all

  for (cs_lnum_t i = 0; i < cs_glob_mesh->n_cells; i++) {
Z[i] = cell_cen[i][2]
        }

Is there anything wrong with the first code or loop?
Thanks

[Yvan Fournier]

Hello,

I assume Z is a field you created. You could simply loop over cells to fill it, as values will be overwritten by the latest contribution. Filling a field containing "Z2-Z1 might be more useful for checking.

There is one very naïve assumption in your code, which the compiler probably warns you about: comparing equality for floating-pont values is not recommendend, as truncation errors will lead to values not being exact.

Instead of "z1 == z2", use something like "fabs(z1 - z2) <epsilon".

Best regards,

Yvan

[Kevin]

Hi. Thanks admin.
That solved it!

Now i wanna ask you two more questions.

1_ Is there any function to get the adjacent or connected boundary face to an interior face?

2_ Is there any function to get the adjacent faces on a boundary face just like cs_glob_mesh_adjacencies which gives the adjacent cells for a cell?

[Yvan Fournier]

Hello,

Are you looking for adjacency with a common edge or with a common cell ? In the second case there might be structures used by CDO whcih are not built by default but could be generated.

Regards,

Yvan

[Kevin]

Hi & thanks.

I need adjacency with common edge.
Sorry, I don't know what is CDO.

[Yvan Fournier]

Hello,

CDO is "Compatible Discrete Operators", an new class of finite volume algorithms being developed in code_saturne.

We have various structures where face adjacencies using edges are built or dedectured susing temporary structures (see in src/mesh/cs_mesh_refine.c for example), but nothing "macro" and simple to use like the permanent structures.

Best regards,

Yvan