Fortran routine: compute the biggest edge of each cell from the biggest edge of each face

[Boone11]

Hi all,

I develop a new turbulence model using Code_Saturne 5.3.4 (the moving to the lastest version will be performed in the near future). I use the fortran routines and I need to compute the longest edge of each cell of my mesh to use it in my turbulence model (the diagonales should not be taken into account). The user function used to do this has to work for hybrid meshes and in the case of non conformal mesh. Furthermore, I would like to avoid creating large table to store data.

I succeed in computing the biggest edge of each face (excluding the diagonales) but I do not know how to link each face to its associated cell to compute the biggest edge of each cell over all its related faces. At the end of the code I developed, I get something like:

Code: Select all

do ifac = 1, nfac		 ! loop over all the faces
	Dmax_face(ifac)=biggest_edge(ifac)		 ! computation of the biggest edge of each face
enddo

I would like to tranform the vector

Code: Select all

Dmax_face(nfac)

into a table that contain the lenght of the edge for all the faces associated to each cell :

Code: Select all

Dmax_face(ncel, nfaces_cell) 		! the faces associated to the current cell (for instance, nfaces_cell=6 if hexa) 

to compute the maximum length associated to each cell as follows:

Code: Select all

Dmax(:)=0.d0
do iel = 1, ncel 		! loop over all the cells
	do i = 1, nfaces_cell 		! loop over of the faces associated to the current cell
		if (Dmax(iel) < Dmax_face(iel,i)) then 		! get the maximum edge length over the faces of the current cell
			Dmax(iel) = Dmax_face(iel,i)
		endif
	enddo
enddo

Is that possible? How can I do this? If not what are the others solutions?
Thanks a lot for your help!

[Yvan Fournier]

Hello,

You do not need extra connectivity, just 3 work arrays, (one on each type of faces, one on cells).

- For each face, loop on edges (defined as the segments joining 2 successive vertices), and compute longest one. In Fortran, use the ipnfac/nodfac indexed arrays for interior faces, ipnfbr/nodfbr for boundary faces. In C, these are the i_face_vtx_idx/i_face_vtx_lst and b_face_vtx_idx/b_face_vtx_lst mameber arrays of the xs_glob_mesh structure.

- Initialize work (results) array on cells to 0.

- For each type of face, use the ifacel array to obtain cell numbers adjacent to an interior face, and ifabor for boundary faces. Update the max length values for the matching cells. In C, these ar i_face_cells and b_face_cells (in the mesh structure).

Regards,

Yvan

[Boone11]

Hi Yvan and thanks a lot for your answer. Here is the code for internal faces, may you just have a look to it and tell me if you see major issues or clumsy loops?

Code: Select all

! FOR INTERNAL FACES

allocate(Lf(nfac)) ! Vector containing the maximum length of each face
allocate(DDmax(ncel)) ! Vector containing the maximum edge length of each cell
allocate(iel1(nfac)) ! Vector containing the index of the node of each face

Lf(:)=0.d0
DDmax(:)=0.d0
iel1(:)=0

!!! Get the maximum edge length of each face and store the node index related to each face
do ifac = 1, nfac ! Loop over the faces
	do ii = ipnfac(ifac),ipnfac(ifac+1)-1 ! Loop to get the position of the nodes of the current internal face
		Vi=nodefac(ii) ! Index of vertice i
		Vip1=nodefac(ii)+1 ! Index of vertice i+1
		! Get the vector coordinates
		S_x=xyznod(1,Vip1)-xyznod(1,Vi)
		S_y=xyznod(2,Vip1)-xyznod(2,Vi)
		S_z=xyznod(3,Vip1)-xyznod(3,Vi)
		L=sqrt(S_x**2+S_y**2+S_z**2) ! Compute the length of the current edge
		if (Lf(ifac) .LT. L) then
			Lf(ifac) = L ! Keep the maximum length edge of each face 
		endif
	enddo
	! Get the cell number associated to the current face 
	iel1(ifac) = ifacel(1,ifac) ! iel1: table containing the index of the node of each face
enddo
!!!

!!! Link the values contained in the table Lf to the related cell and get the maximum edge length of each cell 
do  iel = 1,ncel ! Loop over the cells 
	nface_per_cell=0
	nface_actuelle=1
	do kk=1,nfac
		if (iel == iel1(kk)) then ! If the current cell number match with the cell number indicated by the table iel1
			nface_per_cell=nface_per_cell+1 ! We count the number of face per cell 
		endif 
	enddo 
	do while ((nface_actuelle+nface_per_cell) .LE. nfac)
		do jj=1,nface_per_cell
			if (DDmax(iel) < Lf(jj+nface_actuelle)) then
				DDmax(iel) = Lf(jj+nface_actuelle)  ! We keep the maximum edge length 
			endif 
		enddo 
	enddo 
	nface_actuelle=nface_actuelle+nface_per_cell
enddo 	

Can you confirm me that two successive vertices necessarly make an edge of a cell (as long as the mesh does not contain void)?

Thanks again!

[Yvan Fournier]

Hello,

I did not check your code yet but can confirm 2 succesive vertices defi.e an edge. Be careful also not to forger the last edge, which joins a face's last and first vertex (closing the loop).

Best regards,

Yvan

[Boone11]

Indeed I forgot to close the loop. Thanks for your advices!

I just check Code_Saturne results on Paraview and the the point data ID of the cells vertices. I noticed that two vertices of the same face are not sequential. So the point data ID numerotation in Paraview is different from the cell vertice numerotation in Code_Saturne?

[Yvan Fournier]

Hello,

If you just look at a boundary face mesh for example I would expect the numbering to be similar (though face can be ordered differently, as they are regrouped per face type for postprocessing mesh export).

For cells, things can be different: the cells exported to ParaView use a cell to node connectivity, which is rebuilt by code_saturne, while the internal connectivity is based on face to cells + cells to vertices.

So it really depends on how you test this, but in most cases, ordering will not match.

Best regards,

Yvan

[Boone11]

Ok thanks!

Best regards