/*============================================================================ * Initialization of a nodal connectivity definition based upon * a (possibly partial) descending connectivity *============================================================================*/ /* This file is part of the "Finite Volume Mesh" library, intended to provide finite volume mesh and associated fields I/O and manipulation services. Copyright (C) 2004-2006 EDF This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA */ /*---------------------------------------------------------------------------- * Standard C library headers *----------------------------------------------------------------------------*/ #include #include #include #include /*---------------------------------------------------------------------------- * BFT library headers *----------------------------------------------------------------------------*/ #include #include #include /*---------------------------------------------------------------------------- * Local headers *----------------------------------------------------------------------------*/ #include #include #include #include #include /*---------------------------------------------------------------------------- * Local headers associated with the current file *----------------------------------------------------------------------------*/ #include /*----------------------------------------------------------------------------*/ #ifdef __cplusplus extern "C" { #if 0 } /* Fake brace to force back Emacs auto-indentation back to column 0 */ #endif #endif /* __cplusplus */ /*============================================================================ * Enumeration definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Nodal connectivity reconstruction return code. *----------------------------------------------------------------------------*/ typedef enum { FVM_NODAL_FROM_DESC_SUCCESS, /* Successful reconstruction */ FVM_NODAL_FROM_DESC_FAILURE, /* Reconstruction failure (connectivity pb.) */ FVM_NODAL_FROM_DESC_WORIENT /* Reconstruction with orientation warning (for possible orientation problem) */ } fvm_nodal_from_desc_t; /*============================================================================ * Macro definitions *============================================================================*/ /*============================================================================ * Global static variables *============================================================================*/ /*============================================================================ * Private function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Determine if a given cell is a prism or a polyhedron * * This function should only be called on cells having 2 triangular and * 3 quadrangular faces; if the triangles are adjacent, the cell * is not a prism, and is considered to be a polyhedron. Otherwise, it * seems to be truly a prism (supposing it is closed and well oriented). * * parameters: * cell_id, <-- cell id (0 to n-1) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number index shifts; * size: n_face_lists * face_vertex_idx <-- face -> vertex indexes (per face list) * face_vertex_num <-- face -> vertex numbers (per face list) * cell_face_idx <-- cell -> face indexes (1 to n) * cell_face_num <-- cell -> face numbers (1 to n) * * returns: * type of cell defined by cell_id *----------------------------------------------------------------------------*/ inline static fvm_element_t _is_prism_or_poly(const fvm_lnum_t cell_id, const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[], const fvm_lnum_t *const face_vertex_num[], const fvm_lnum_t cell_face_idx[], const fvm_lnum_t cell_face_num[]) { int vtx_id_1, vtx_id_2; fvm_lnum_t face_id, fl; fvm_lnum_t idx, idx_start, idx_end; fvm_lnum_t n_face_vertices; fvm_lnum_t vtx, vertex_id_start, vertex_id_end; fvm_lnum_t vtx_tria[6]; int n_trias = 0; /* Extract 2 triangles */ idx_start = cell_face_idx[cell_id] - 1; idx_end = cell_face_idx[cell_id + 1] - 1; for (idx = idx_start ; idx < idx_end ; idx++) { face_id = FVM_ABS(cell_face_num[idx]) - 1; for (fl = n_face_lists - 1 ; face_id < face_list_shift[fl] ; fl--); assert(fl > -1); face_id -= face_list_shift[fl]; vertex_id_start = face_vertex_idx[fl][face_id] - 1; vertex_id_end = face_vertex_idx[fl][face_id + 1] - 1; n_face_vertices = vertex_id_end - vertex_id_start; if (n_face_vertices == 3) { if (cell_face_num[idx] > 0) { for (vtx = 0 ; vtx < 3 ; vtx++) vtx_tria[n_trias*3 + vtx] = face_vertex_num[fl][vertex_id_start + vtx]; } else { for (vtx = 0 ; vtx < 3 ; vtx++) vtx_tria[n_trias*3 + vtx] = face_vertex_num[fl][vertex_id_end - 1 - vtx]; } n_trias += 1; if (n_trias == 2) /* Once we have found the two triangles, */ break; /* we do not need to look at other faces */ } } assert(n_trias == 2); /* Are the triangles adjacent ? */ for (vtx_id_1 = 0; vtx_id_1 < 3; vtx_id_1++) { for (vtx_id_2 = 0; vtx_id_2 < 3; vtx_id_2++) { if (vtx_tria[3 + vtx_id_2] == vtx_tria[vtx_id_1]) { return FVM_CELL_POLY; } } } /* If no adjacent triangles were found, we have a prism */ return FVM_CELL_PRISM; } /*---------------------------------------------------------------------------- * Determination of a given cell's type. * * If the optional cell_vtx_tria[3*4] and cell_vtx_quad[4*6] arrays are given, * they are filled with the vertex indexes of the cell's triangle and * quadrangle type faces (as long as no simple polyhedral face is encountered * and we have no more than 4 triangles or 6 quadrangles), for easier * nodal cell reconstuction. Vertex indexes are ordered as usual for * an outward pointing normal. * * parameters: * cell_id, <-- cell id (0 to n-1) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number index shifts; * size: n_face_lists * face_vertex_idx <-- face -> vertex indexes (per face list) * face_vertex_num <-- face -> vertex numbers (per face list) * cell_face_idx <-- cell -> face indexes (1 to n) * cell_face_num <-- cell -> face numbers (1 to n) * cell_vtx_tria --> local triangle definitions (optional, 4 max.) * cell_vtx_quad --> local quadrangle definitions (optional, 6 max.) * * returns: * type of cell defined by cell_id *----------------------------------------------------------------------------*/ inline static fvm_element_t _nodal_cell_from_desc(const fvm_lnum_t cell_id, const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[], const fvm_lnum_t *const face_vertex_num[], const fvm_lnum_t cell_face_idx[], const fvm_lnum_t cell_face_num[], fvm_lnum_t *const cell_vtx_tria, fvm_lnum_t *const cell_vtx_quad) { fvm_lnum_t n_cell_faces; fvm_lnum_t face_id, fl; fvm_lnum_t idx, idx_start, idx_end; fvm_lnum_t n_face_vertices; fvm_lnum_t vtx, vertex_id_start, vertex_id_end; fvm_element_t cell_type; fvm_lnum_t n_trias = 0; fvm_lnum_t n_quads = 0; fvm_lnum_t n_ngons = 0; /* Guess connectivity types */ /*--------------------------*/ n_cell_faces = cell_face_idx[cell_id + 1] - cell_face_idx[cell_id]; /* If we have more than 6 faces, we have a general simple polyhedron */ if (n_cell_faces > 6) return FVM_CELL_POLY; /* Otherwise, we probably have a "classical" element; for example, with 6 faces, we probably have a hexahedron, though we could also have a tetrahedron with one of its edges and thus two of its faces split. We count vertices per face to check. */ idx_start = cell_face_idx[cell_id] - 1; idx_end = cell_face_idx[cell_id + 1] - 1; for (idx = idx_start ; idx < idx_end ; idx++) { face_id = FVM_ABS(cell_face_num[idx]) - 1; for (fl = n_face_lists - 1 ; face_id < face_list_shift[fl] ; fl--); assert(fl > -1); face_id -= face_list_shift[fl]; vertex_id_start = face_vertex_idx[fl][face_id] - 1; vertex_id_end = face_vertex_idx[fl][face_id + 1] - 1; n_face_vertices = vertex_id_end - vertex_id_start; if (n_face_vertices == 3) { if (cell_vtx_tria != NULL && n_trias < 4) { if (cell_face_num[idx] > 0) { for (vtx = 0 ; vtx < n_face_vertices ; vtx++) cell_vtx_tria[n_trias*3 + vtx] = face_vertex_num[fl][vertex_id_start + vtx]; } else { for (vtx = 0 ; vtx < n_face_vertices ; vtx++) cell_vtx_tria[n_trias*3 + vtx] = face_vertex_num[fl][vertex_id_end - 1 - vtx]; } } n_trias += 1; } else if (n_face_vertices == 4) { if (cell_vtx_quad != NULL && n_quads < 6) { if (cell_face_num[idx] > 0) { for (vtx = 0 ; vtx < n_face_vertices ; vtx++) cell_vtx_quad[n_quads*4 + vtx] = face_vertex_num[fl][vertex_id_start + vtx]; } else { for (vtx = 0 ; vtx < n_face_vertices ; vtx++) cell_vtx_quad[n_quads*4 + vtx] = face_vertex_num[fl][vertex_id_end - 1 - vtx]; } } n_quads += 1; } else n_ngons += 1; } /* Return element type */ if (n_ngons > 0) cell_type = FVM_CELL_POLY; else { if (n_trias == 0 && n_quads == 6) cell_type = FVM_CELL_HEXA; else if (n_trias == 2 && n_quads == 3) cell_type = _is_prism_or_poly(cell_id, n_face_lists, face_list_shift, face_vertex_idx, face_vertex_num, cell_face_idx, cell_face_num); else if (n_trias == 4) { if (n_quads == 0) cell_type = FVM_CELL_TETRA; else if (n_quads == 1) cell_type = FVM_CELL_PYRAM; else cell_type = FVM_CELL_POLY; } else cell_type = FVM_CELL_POLY; } return cell_type; } /*---------------------------------------------------------------------------- * Tetrahedron construction * * parameters: * cell_vtx_tria <-- triangular faces connectivity * cell_vtx_tetra --> tetrahedron connectivity (pre-allocated) *----------------------------------------------------------------------------*/ inline static fvm_nodal_from_desc_t _nodal_from_desc_cnv_cel_tetra(const fvm_lnum_t cell_vtx_tria[], fvm_lnum_t cell_vtx_tetra[]) { fvm_lnum_t vertex_id, face_id; fvm_lnum_t direction; fvm_lnum_t vtx_num, vtx_num_1, vtx_num_2; _Bool warn_orient = false; #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("face 1 : %d %d %d\n", cell_vtx_tria[0], cell_vtx_tria[1], cell_vtx_tria[2]); bft_printf("face 2 : %d %d %d\n", cell_vtx_tria[3], cell_vtx_tria[4], cell_vtx_tria[5]); bft_printf("face 3 : %d %d %d\n", cell_vtx_tria[6], cell_vtx_tria[7], cell_vtx_tria[8]); bft_printf("face 4 : %d %d %d\n", cell_vtx_tria[9], cell_vtx_tria[10], cell_vtx_tria[11]); #endif /* Base : vertices of the first face; we take the vertices in opposite order ("bottom" face numbering of the tetrahedron with outward pointing normal). * x 4 * /|\ * / | \ * / | \ * 1 x- -|- -x 3 * \ | / * \ | / * \|/ * x 2 */ cell_vtx_tetra[0] = cell_vtx_tria[2]; cell_vtx_tetra[1] = cell_vtx_tria[1]; cell_vtx_tetra[2] = cell_vtx_tria[0]; /* We have found 3 of 4 vertices; all other triangles should share vertex 4, and one of those should share vertices 1 and 2 of the base triangle. */ vtx_num_1 = cell_vtx_tetra[0]; vtx_num_2 = cell_vtx_tetra[1]; direction = 0; for (face_id = 1 ; face_id < 4 ; face_id++) { for (vertex_id = 0 ; vertex_id < 3 ; vertex_id++) { vtx_num = cell_vtx_tria[face_id*3 + vertex_id]; if (vtx_num == vtx_num_1) { if (cell_vtx_tria[face_id*3 + ((vertex_id+1) % 3)] == vtx_num_2) { direction = 1; break; } else if (cell_vtx_tria[face_id*3 + ((vertex_id-1+3) % 3)] == vtx_num_2) { direction = -1; break; } } } if (direction != 0) break; } if (direction == -1) warn_orient = true; else if (direction == 0) return FVM_NODAL_FROM_DESC_FAILURE; cell_vtx_tetra[3] = cell_vtx_tria[face_id*3 + ((vertex_id + (3 + (2 * direction))) % 3)]; #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("tetra : %d %d %d %d\n", cell_vtx_tetra[0], cell_vtx_tetra[1], cell_vtx_tetra[2], cell_vtx_tetra[3]); #endif if (warn_orient == true) return FVM_NODAL_FROM_DESC_WORIENT; else return FVM_NODAL_FROM_DESC_SUCCESS; } /*---------------------------------------------------------------------------- * Pyramid construction * * parameters: * cell_vtx_tria <-- triangular faces connectivity * cell_vtx_quad <-- quadrangle faces connectivity * cell_vtx_pyram --> pyramid connectivity (pre-allocated) *----------------------------------------------------------------------------*/ inline static fvm_nodal_from_desc_t _nodal_from_desc_cnv_cel_pyram(const fvm_lnum_t cell_vtx_tria[], const fvm_lnum_t cell_vtx_quad[], fvm_lnum_t cell_vtx_pyram[]) { fvm_lnum_t vertex_id, face_id; fvm_lnum_t direction; fvm_lnum_t vtx_num, vtx_num_1, vtx_num_2; _Bool warn_orient = false; #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("face 1 : %d %d %d %d\n", cell_vtx_quad[0], cell_vtx_quad[1], cell_vtx_quad[2], cell_vtx_quad[3]); bft_printf("face 2 : %d %d %d\n", cell_vtx_tria[0], cell_vtx_tria[1], cell_vtx_tria[2]); bft_printf("face 3 : %d %d %d\n", cell_vtx_tria[3], cell_vtx_tria[4], cell_vtx_tria[5]); bft_printf("face 4 : %d %d %d\n", cell_vtx_tria[6], cell_vtx_tria[7], cell_vtx_tria[8]); bft_printf("face 5 : %d %d %d\n", cell_vtx_tria[9], cell_vtx_tria[10], cell_vtx_tria[11]); #endif /* Base : vertices of the quadrangle; we take the vertices in opposite order ("bottom" face numbering of the pyramid with outward pointing normal). * 5 x * /|\ * //| \ * // | \ * 4 x/--|---x 3 * // | / * // | / * 1 x-------x 2 */ cell_vtx_pyram[0] = cell_vtx_quad[3]; cell_vtx_pyram[1] = cell_vtx_quad[2]; cell_vtx_pyram[2] = cell_vtx_quad[1]; cell_vtx_pyram[3] = cell_vtx_quad[0]; /* We have found 4 out of 5 vertices; all 4 triangles should share vertex 5, and one of those should share vertices 1 and 2 of the base quadrangle. */ vtx_num_1 = cell_vtx_pyram[0]; vtx_num_2 = cell_vtx_pyram[1]; direction = 0; for (face_id = 0 ; face_id < 4 ; face_id++) { for (vertex_id = 0 ; vertex_id < 3 ; vertex_id++) { vtx_num = cell_vtx_tria[face_id*3 + vertex_id]; if (vtx_num == vtx_num_1) { if (cell_vtx_tria[face_id*3 + ((vertex_id+1) % 3)] == vtx_num_2) { direction = 1; break; } else if (cell_vtx_tria[face_id*3 + ((vertex_id-1+3) % 3)] == vtx_num_2) { direction = -1; break; } } } if (direction != 0) break; } if (direction == -1) warn_orient = true; else if (direction == 0) return FVM_NODAL_FROM_DESC_FAILURE; cell_vtx_pyram[4] = cell_vtx_tria[face_id*3 + ((vertex_id + (3 + (2 * direction))) % 3)]; #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("pyram : %d %d %d %d %d\n", cell_vtx_pyram[0], cell_vtx_pyram[1], cell_vtx_pyram[2], cell_vtx_pyram[3], cell_vtx_pyram[4]); #endif if (warn_orient == true) return FVM_NODAL_FROM_DESC_WORIENT; else return FVM_NODAL_FROM_DESC_SUCCESS; } /*---------------------------------------------------------------------------- * Prism (pentahedron) construction * * parameters: * cell_vtx_tria <-- triangular faces connectivity * cell_vtx_quad <-- quadrangle faces connectivity * cell_vtx_prism --> prism connectivity (pre-allocated) *----------------------------------------------------------------------------*/ inline static fvm_nodal_from_desc_t _nodal_from_desc_cnv_cel_prism(const fvm_lnum_t cell_vtx_tria[], const fvm_lnum_t cell_vtx_quad[], fvm_lnum_t cell_vtx_prism[]) { fvm_lnum_t vertex_id, face_id; fvm_lnum_t ipass, direction; fvm_lnum_t vtx_num, vtx_num_1, vtx_num_2; _Bool warn_orient = false; #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("face 1 : %d %d %d\n", cell_vtx_tria[0], cell_vtx_tria[1], cell_vtx_tria[2]); bft_printf("face 2 : %d %d %d\n", cell_vtx_tria[3], cell_vtx_tria[4], cell_vtx_tria[5]); bft_printf("face 3 : %d %d %d %d\n", cell_vtx_quad[0], cell_vtx_quad[1], cell_vtx_quad[2], cell_vtx_quad[3]); bft_printf("face 4 : %d %d %d %d\n", cell_vtx_quad[4], cell_vtx_quad[5], cell_vtx_quad[6], cell_vtx_quad[7]); bft_printf("face 5 : %d %d %d %d\n", cell_vtx_quad[8], cell_vtx_quad[9], cell_vtx_quad[10], cell_vtx_quad[11]); #endif /* Base : vertices of the first triangle; we take the vertices in opposite order ("bottom" face numbering of the prism with outward pointing normal). * 4 x-------x 6 * |\ /| * | \ / | * 1 x- \-/ -x 3 * \ 5x / * \ | / * \|/ * x 2 */ cell_vtx_prism[0] = cell_vtx_tria[2]; cell_vtx_prism[1] = cell_vtx_tria[1]; cell_vtx_prism[2] = cell_vtx_tria[0]; /* We have found 3 out of 6 vertices; we first seek the quadrangle sharing vertices 1 and 2, so as to determine vertices 4 and 5, then we seek the quadrangle sharing vertices 2 and 3, so as to determine vertices 5 and 6. */ for (ipass = 0 ; ipass < 2 ; ipass++) { vtx_num_1 = cell_vtx_prism[ ipass]; vtx_num_2 = cell_vtx_prism[1 + ipass]; direction = 0; for (face_id = 0 ; face_id < 4 ; face_id++) { for (vertex_id = 0 ; vertex_id < 4 ; vertex_id++) { vtx_num = cell_vtx_quad[face_id*4 + vertex_id]; if (vtx_num == vtx_num_1) { if (cell_vtx_quad[face_id*4 + ((vertex_id+1) % 4)] == vtx_num_2) { direction = 1; break; } else if ( cell_vtx_quad[face_id*4 + ((vertex_id-1+4) % 4)] == vtx_num_2) { direction = -1; break; } } } if (direction != 0) break; } if (direction == -1) warn_orient = true; else if (direction == 0) return FVM_NODAL_FROM_DESC_FAILURE; cell_vtx_prism[3 + ipass] = cell_vtx_quad[face_id*4 + ((vertex_id + (4 + (3 * direction))) % 4)]; cell_vtx_prism[4 + ipass] = cell_vtx_quad[face_id*4 + ((vertex_id + (4 + (2 * direction))) % 4)]; } #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("prism : %d %d %d %d %d %d\n", cell_vtx_prism[0], cell_vtx_prism[1], cell_vtx_prism[2], cell_vtx_prism[3], cell_vtx_prism[4], cell_vtx_prism[5]); #endif if (warn_orient == true) return FVM_NODAL_FROM_DESC_WORIENT; else return FVM_NODAL_FROM_DESC_SUCCESS; } /*---------------------------------------------------------------------------- * Hexahedron construction * * parameters: * cell_vtx_quad <-- quadrangle faces connectivity * cell_vtx_hexa --> hexahedron connectivity (pre-allocated) *----------------------------------------------------------------------------*/ inline static fvm_nodal_from_desc_t _nodal_from_desc_cnv_cel_hexa(const fvm_lnum_t cell_vtx_quad[], fvm_lnum_t cell_vtx_hexa[]) { fvm_lnum_t vertex_id, face_id; fvm_lnum_t ipass, direction; fvm_lnum_t vtx_num, vtx_num_1, vtx_num_2; _Bool warn_orient = false; #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("face 1 : %d %d %d %d\n", cell_vtx_quad[0], cell_vtx_quad[1], cell_vtx_quad[2], cell_vtx_quad[3]); bft_printf("face 2 : %d %d %d %d\n", cell_vtx_quad[4], cell_vtx_quad[5], cell_vtx_quad[6], cell_vtx_quad[7]); bft_printf("face 3 : %d %d %d %d\n", cell_vtx_quad[8], cell_vtx_quad[9], cell_vtx_quad[10], cell_vtx_quad[11]); bft_printf("face 4 : %d %d %d %d\n", cell_vtx_quad[12], cell_vtx_quad[13], cell_vtx_quad[14], cell_vtx_quad[15]); bft_printf("face 5 : %d %d %d %d\n", cell_vtx_quad[16], cell_vtx_quad[17], cell_vtx_quad[18], cell_vtx_quad[19]); bft_printf("face 6 : %d %d %d %d\n", cell_vtx_quad[20], cell_vtx_quad[21], cell_vtx_quad[22], cell_vtx_quad[23]); #endif /* Base : vertices of the fisrt face; we take the vertices in opposite order ("bottom" face numbering of the pyramid with outward pointing normal). * 8 x-------x 7 * /| /| * / | / | * 5 x-------x6 | * | 4x----|--x 3 * | / | / * |/ |/ * 1 x-------x 2 */ cell_vtx_hexa[0] = cell_vtx_quad[3]; cell_vtx_hexa[1] = cell_vtx_quad[2]; cell_vtx_hexa[2] = cell_vtx_quad[1]; cell_vtx_hexa[3] = cell_vtx_quad[0]; /* We have found 4 out of 8 vertices; we first seek the quadrangle sharing vertices 1 and 2, so as to determine vertices 5 and 6, then we seek the quadrangle sharing vertices 3 and 4, so as to determine vertices 7 and 8. */ for (ipass = 0 ; ipass < 2 ; ipass++) { vtx_num_1 = cell_vtx_hexa[ ipass * 2 ]; vtx_num_2 = cell_vtx_hexa[1 + (ipass * 2)]; direction = 0; for (face_id = 1 ; face_id < 6 ; face_id++) { for (vertex_id = 0 ; vertex_id < 4 ; vertex_id++) { vtx_num = cell_vtx_quad[face_id*4 + vertex_id]; if (vtx_num == vtx_num_1) { if (cell_vtx_quad[face_id*4 + ((vertex_id+1) % 4)] == vtx_num_2) { direction = 1; break; } else if ( cell_vtx_quad[face_id*4 + ((vertex_id-1+4) % 4)] == vtx_num_2) { direction = -1; break; } } } if (direction != 0) break; } if (direction == -1) warn_orient = true; else if (direction == 0) return FVM_NODAL_FROM_DESC_FAILURE; cell_vtx_hexa[4 + (ipass * 2)] = cell_vtx_quad[face_id*4 + ((vertex_id + (4 + (3 * direction))) % 4)]; cell_vtx_hexa[5 + (ipass * 2)] = cell_vtx_quad[face_id*4 + ((vertex_id + (4 + (2 * direction))) % 4)]; } #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("hexa : %d %d %d %d %d %d %d %d\n", cell_vtx_hexa[0], cell_vtx_hexa[1], cell_vtx_hexa[2], cell_vtx_hexa[3], cell_vtx_hexa[4], cell_vtx_hexa[5], cell_vtx_hexa[6], cell_vtx_hexa[7]); #endif if (warn_orient == true) return FVM_NODAL_FROM_DESC_WORIENT; else return FVM_NODAL_FROM_DESC_SUCCESS; } /*---------------------------------------------------------------------------- * Determination of the number of vertices defining a given face. * * parameters: * face_id, <-- face id (0 to n-1) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number index shifts; * size: n_face_lists * face_vertex_idx <-- face -> vertex indexes (per face list) * * returns: * number of vertices of face defined by face_id *----------------------------------------------------------------------------*/ inline static fvm_lnum_t _nodal_face_from_desc_size(const fvm_lnum_t face_id, const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[]) { fvm_lnum_t fl, _face_id; fvm_lnum_t n_face_vertices; fvm_lnum_t vertex_id_start, vertex_id_end; /* Compute number of vertices */ /*----------------------------*/ _face_id = face_id; for (fl = n_face_lists - 1 ; _face_id < face_list_shift[fl] ; fl--); assert(fl > -1); _face_id -= face_list_shift[fl]; vertex_id_start = face_vertex_idx[fl][_face_id] - 1; vertex_id_end = face_vertex_idx[fl][_face_id + 1] - 1; n_face_vertices = vertex_id_end - vertex_id_start; return n_face_vertices; } /*---------------------------------------------------------------------------- * Copy the vertex numbers defining a given face. * * The face_vtx[] array is filled with the vertex numbers of the face, * and should be large enough to receive them. * * parameters: * face_id, <-- face id (0 to n-1) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number index shifts; * size: n_face_lists * face_vertex_idx <-- face -> vertex indexes (per face list) * face_vertex_num <-- face -> vertex numbers (per face list) * face_vtx --> local face definition *----------------------------------------------------------------------------*/ inline static void _nodal_face_from_desc_copy(const fvm_lnum_t face_id, const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[], const fvm_lnum_t *const face_vertex_num[], fvm_lnum_t *const face_vtx) { fvm_lnum_t fl, _face_id; fvm_lnum_t vtx, vertex_id, vertex_id_start, vertex_id_end; /* Copy vertex numbers */ /*---------------------*/ _face_id = face_id; for (fl = n_face_lists - 1 ; _face_id < face_list_shift[fl] ; fl--); assert(fl > -1); _face_id -= face_list_shift[fl]; vertex_id_start = face_vertex_idx[fl][_face_id] - 1; vertex_id_end = face_vertex_idx[fl][_face_id + 1] - 1; for (vtx = 0, vertex_id = vertex_id_start ; vertex_id < vertex_id_end ; vtx++, vertex_id++) face_vtx[vtx] = face_vertex_num[fl][vertex_id]; } /*---------------------------------------------------------------------------- * Polyhedral cells connectivity extraction. * * parameters: * this_nodal <-> nodal mesh structure * n_polys <-- size of list_poly[] * list_poly <-- list of polyhedral cells (cell index, 1 to n) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number number index shifts; * size: n_face_lists * face_vertex_idx <-- face -> vertex indexes (per face list) * face_vertex_num <-- face -> vertex numbers (per face list) * cell_face_idx <-- cell -> face indexes (1 to n) * cell_face_num <-- cell -> face numbers (1 to n) * cell_face_list --> numbers of faces defining polyhedra *----------------------------------------------------------------------------*/ static void _fvm_nodal_extract_polyhedra(fvm_nodal_section_t *const this_section, const fvm_lnum_t n_polys, const fvm_lnum_t list_poly[], const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[], const fvm_lnum_t *const face_vertex_num[], const fvm_lnum_t cell_face_idx[], const fvm_lnum_t cell_face_num[], fvm_lnum_t * cell_face_list[]) { fvm_lnum_t n_faces, n_cell_faces; fvm_lnum_t c_cell_face_vertex_idxs, c_cell_face_vertex_nums; fvm_lnum_t n_cell_face_vertex_nums; fvm_lnum_t face_counter, face_id, poly_id, cell_id; fvm_lnum_t fl, i, idx, idx_start, idx_end, num_count; fvm_lnum_t *local_face_num = NULL; int sgn; /* Indicators initialization */ n_faces = face_list_shift[n_face_lists] - face_list_shift[0]; BFT_MALLOC(local_face_num, n_faces, fvm_lnum_t); for (face_id = 0 ; face_id < n_faces ; face_id++) local_face_num[face_id] = 0; /* Flagging of referenced faces and Cells -> Faces indexes */ /*---------------------------------------------------------*/ n_cell_faces = 0; BFT_MALLOC(this_section->_face_index, n_polys + 1, fvm_lnum_t); this_section->face_index = this_section->_face_index; this_section->_face_index[0] = 0; for (poly_id = 0 ; poly_id < n_polys ; poly_id++) { cell_id = list_poly[poly_id] - 1; idx_start = cell_face_idx[cell_id] - 1; idx_end = cell_face_idx[cell_id +1] - 1; this_section->_face_index[poly_id + 1] = this_section->_face_index[poly_id] + (idx_end - idx_start); for (idx = idx_start ; idx < idx_end ; idx++) { face_id = FVM_ABS(cell_face_num[idx]) - 1; /* Mark only used values for now, local_face_num[] values will be computed later by looping on faces so that for nonzero values, local_face_num[i] > local_face_num[j] if i > j ; this is important for later parts of this algorithm */ if (local_face_num[face_id] == 0) local_face_num[face_id] = 1; } } /* Counting for faces -> vertices connectivity and local face numbering */ n_cell_face_vertex_nums = 0; for (face_counter = 0 ; face_counter < n_faces ; face_counter++) { if (local_face_num[face_counter] != 0) { /* Transform local_face_num[] from a marker to a renumbering array */ n_cell_faces += 1; local_face_num[face_counter] = n_cell_faces; /* Counting for faces -> vertices connectivity */ face_id = face_counter; for (fl = n_face_lists - 1 ; face_id < face_list_shift[fl] ; fl--); assert(fl > -1); face_id -= face_list_shift[fl]; n_cell_face_vertex_nums += ( face_vertex_idx[fl][face_id+1] - face_vertex_idx[fl][face_id]); } } /* Cells -> Faces Connectivity (face numbers) */ /*--------------------------------------------*/ BFT_MALLOC(this_section->_face_num, this_section->_face_index[n_polys], fvm_lnum_t); this_section->face_num = this_section->_face_num; num_count = 0; for (poly_id = 0 ; poly_id < n_polys ; poly_id++) { cell_id = list_poly[poly_id] - 1; idx_start = cell_face_idx[cell_id] - 1; idx_end = cell_face_idx[cell_id +1] - 1; for (idx = idx_start ; idx < idx_end ; idx++) { face_id = FVM_ABS(cell_face_num[idx]) - 1; sgn = (cell_face_num[idx] > 0) ? 1 : -1; this_section->_face_num[num_count++] = sgn * local_face_num[face_id] ; } } assert(num_count == this_section->_face_index[n_polys]); /* Faces -> Vertices Connectivity */ /*--------------------------------*/ if (cell_face_list != NULL) BFT_MALLOC(*cell_face_list, n_cell_faces, fvm_lnum_t); BFT_MALLOC(this_section->_vertex_index, n_cell_faces + 1, fvm_lnum_t); this_section->vertex_index = this_section->_vertex_index; BFT_MALLOC(this_section->_vertex_num, n_cell_face_vertex_nums, fvm_lnum_t); this_section->vertex_num = this_section->_vertex_num; /* Definition */ c_cell_face_vertex_idxs = 0; c_cell_face_vertex_nums = 0; this_section->_vertex_index[0] = 0; for (face_counter = 0 ; face_counter < n_faces ; face_counter++) { if (local_face_num[face_counter] != 0) { if (cell_face_list != NULL) (*cell_face_list)[local_face_num[face_counter] -1 ] = face_counter + 1; face_id = face_counter; for (fl = n_face_lists - 1 ; face_id < face_list_shift[fl] ; fl--); assert(fl < n_face_lists); face_id -= face_list_shift[fl]; for (i = face_vertex_idx[fl][face_id] - 1 ; i < face_vertex_idx[fl][face_id + 1] - 1 ; i++) this_section->_vertex_num[c_cell_face_vertex_nums++] = face_vertex_num[fl][i]; c_cell_face_vertex_idxs += 1; this_section->_vertex_index[c_cell_face_vertex_idxs] = c_cell_face_vertex_nums; } } /* No past-the-end index value counted, so counter = this_nodal->n_cell_faces */ assert(c_cell_face_vertex_idxs == n_cell_faces); assert(c_cell_face_vertex_nums == n_cell_face_vertex_nums); /* Set pointers and connectivity size */ this_section->connectivity_size = n_cell_face_vertex_nums; this_section->n_faces = n_cell_faces; this_section->face_index = this_section->_face_index; this_section->face_num = this_section->_face_num; this_section->vertex_index = this_section->_vertex_index; this_section->vertex_num = this_section->_vertex_num; /* Free memory */ BFT_FREE(local_face_num); } /*---------------------------------------------------------------------------- * Raise parent element numbering in given sections by one level, as * defined by the parent_element_num[] arrays. * * This is useful when a nodal mesh is extracted from a temporary subset of * a parent mesh (corresponding to the parent_element_num[] list, using * 1 to n numbering), and the final nodal mesh element parent numbering * should correspond to that parent mesh and not the temporary subset. * * parameters: * n_sections <-- size of sections array * sections <-- array of sections to reduce * parent_element_num <-- element -> parent element number (1 to n) if * non-trivial (i.e. if element definitions * correspond to a subset of the parent mesh), * NULL otherwise. *----------------------------------------------------------------------------*/ static void _raise_sections_parent_num(const int n_sections, fvm_nodal_section_t *const sections[], const fvm_lnum_t parent_element_num[]) { int section_id; fvm_lnum_t element_counter; fvm_nodal_section_t *section; if (parent_element_num == NULL) return; for (section_id = 0 ; section_id < n_sections ; section_id++) { section = sections[section_id]; if (section != NULL) { if (section->_parent_element_num == NULL) { BFT_MALLOC(section->_parent_element_num, section->n_elements, fvm_lnum_t); section->parent_element_num = section->_parent_element_num; } for (element_counter = 0 ; element_counter < section->n_elements ; element_counter++) section->_parent_element_num[element_counter] = parent_element_num[section->parent_element_num[element_counter] - 1]; } } } /*---------------------------------------------------------------------------- * Free parent element numbering arrays when not needed. * * If the parent mesh is already complete and ordered for a given section * (i.e. not based on a partial extraction or then based on the n first * elements in the same order), the parent element number is not needed, * and may be freed. * * parameters: * sections <-- array of sections to reduce * n_sections <-- size of sections array *----------------------------------------------------------------------------*/ static void _optimize_sections_parent_num(const int n_sections, fvm_nodal_section_t *const sections[]) { int section_id; fvm_lnum_t element_counter; fvm_nodal_section_t *section; /* If the parent mesh is already complete for a given element type (i.e. not based on a partial extraction or then based on the n first elements), the parent cell number is not needed */ for (section_id = 0 ; section_id < n_sections ; section_id++) { section = sections[section_id]; if (section != NULL) { for (element_counter = 0 ; element_counter < section->n_elements ; element_counter++) { if (section->parent_element_num[element_counter] != element_counter + 1) break; } if (element_counter == section->n_elements) { section->parent_element_num = NULL; BFT_FREE(section->_parent_element_num); } } } } /*---------------------------------------------------------------------------- * Add nodal mesh structure sections to a nodal mesh. * * Sections to add are defined by an array of section pointers, * which may contain NULL entries. Only the non-NULL entries are * added to the nodal mesh structure, and they belong to that structure * from this point on (i.e. their lifecycle is based upon that of * the nodal mesh structure). * * parameters: * this_nodal <-> nodal mesh structure * n_sections <-- size of sections array * sections <-- array of sections to add *----------------------------------------------------------------------------*/ static void _fvm_nodal_add_sections(fvm_nodal_t *const this_nodal, const int n_sections, fvm_nodal_section_t *const sections[]) { int section_id, section_count; fvm_nodal_section_t *section; /* Add sections to nodal mesh structure */ section_count = 0; for (section_id = 0 ; section_id < n_sections ; section_id++) { section = sections[section_id]; if (section != NULL) section_count++; } BFT_REALLOC(this_nodal->sections, this_nodal->n_sections + section_count, fvm_nodal_section_t *); section_count = 0; for (section_id = 0 ; section_id < n_sections ; section_id++) { section = sections[section_id]; if (section != NULL) this_nodal->sections[this_nodal->n_sections + section_count++] = section; } this_nodal->n_sections += section_count; } /*============================================================================ * Public function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Convert and add cells from an descending connectivity mesh to a nodal mesh. * * If the optional filter list extr_cells[] argument is non-NULL, cells * {extr_cells[0], extr_cells[1], extr_cells[n_extr_cells - 1]} are converted * and added to the nodal mesh. If this filter is set to NULL, cells * {1, 2, ..., n_extr_cells} are considered. * * In addition, an optional parent_cell_num[] array may also be given, in * case the descending connectivity mesh definition is based on a temporary * subset of a parent mesh, (corresponding to the parent_cell_num[] list, * using 1 to n numbering), and the final nodal mesh element parent numbering * should correspond to that parent mesh and not the temporary subset. * * parameters: * this_nodal <-> nodal mesh structure * n_extr_cells <-- count of cells to add * extr_cells <-- optional filter list of cells to extract (1 to n) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number index shifts; * size: n_face_lists + 1 * face_vertex_idx <-- face -> vertex indexes (per face list) * face_vertex_num <-- face -> vertex numbers (per face list) * cell_face_idx <-- cell -> face indexes (1 to n) * cell_face_num <-- cell -> face numbers (1 to n) * parent_cell_num <-- cell -> parent cell number (1 to n) if non-trivial * (i.e. if cell definitions correspond to a subset * of the parent mesh), NULL otherwise. * cell_face_list --> numbers of faces defining polyhedra *----------------------------------------------------------------------------*/ void fvm_nodal_from_desc_add_cells(fvm_nodal_t *const this_nodal, const fvm_lnum_t n_extr_cells, const fvm_lnum_t extr_cells[], const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[], const fvm_lnum_t *const face_vertex_num[], const fvm_lnum_t cell_face_idx[], const fvm_lnum_t cell_face_num[], const fvm_lnum_t parent_cell_num[], fvm_lnum_t * cell_face_list[]) { int type_id; fvm_lnum_t cell_counter, cell_id; fvm_element_t cell_type; fvm_lnum_t cell_vtx_tria[3*4]; /* We will only seek to fill these arrays */ fvm_lnum_t cell_vtx_quad[4*6]; /* for local faces 1-4 and 1-6 at most */ fvm_lnum_t *p_cell_vertex; fvm_lnum_t n_elements_type[FVM_N_ELEMENT_TYPES]; fvm_gnum_t n_g_elements_type[FVM_N_ELEMENT_TYPES]; fvm_nodal_section_t *section; fvm_nodal_section_t *sections[FVM_N_ELEMENT_TYPES]; fvm_gnum_t n_orient_pbs = 0; /* Number of cells with potential (non- fatal) orientation problem */ fvm_nodal_from_desc_t retcode; /* Initialization */ for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) { n_elements_type[type_id] = 0; sections[type_id] = NULL; } /* Guess connectivity types */ /*--------------------------*/ for (cell_counter = 0 ; cell_counter < n_extr_cells ; cell_counter++) { if (extr_cells != NULL) cell_id = extr_cells[cell_counter] - 1; else cell_id = cell_counter; cell_type = _nodal_cell_from_desc(cell_id, n_face_lists, face_list_shift, face_vertex_idx, face_vertex_num, cell_face_idx, cell_face_num, NULL, NULL); n_elements_type[cell_type] += 1; } #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("dbg : fvm_nodal_cells_from_desc\n" "n_tetra = %d\n" "n_pyram = %d\n" "n_prism = %d\n" "n_hexa = %d\n" "n_poly = %d\n", n_elements_type[FVM_CELL_TETRA], n_elements_type[FVM_CELL_PYRAM], n_elements_type[FVM_CELL_PRISM], n_elements_type[FVM_CELL_HEXA], n_elements_type[FVM_CELL_POLY]); #endif /* Set dimensions (and reset local counters) */ for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) n_g_elements_type[type_id] = n_elements_type[type_id]; fvm_parall_counter(n_g_elements_type, FVM_N_ELEMENT_TYPES); for (cell_type = FVM_CELL_TETRA ; cell_type <= FVM_CELL_POLY ; cell_type++) { if (n_g_elements_type[cell_type] > 0) { sections[cell_type] = fvm_nodal_section_create(cell_type); sections[cell_type]->n_elements = n_elements_type[cell_type]; this_nodal->n_cells += n_elements_type[cell_type]; } n_elements_type[cell_type] = 0; } /* Main memory allocations */ for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) { section = sections[type_id]; if (section != NULL) { if ( section->type != FVM_FACE_POLY && section->type != FVM_CELL_POLY) { section->stride = fvm_nodal_n_vertices_element[type_id]; section->connectivity_size = section->stride * section->n_elements; BFT_MALLOC(section->_vertex_num, section->connectivity_size, fvm_lnum_t); section->vertex_num = section->_vertex_num; } } } for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) { section = sections[type_id]; if (section != NULL) { BFT_MALLOC(section->_parent_element_num, section->n_elements, fvm_lnum_t); section->parent_element_num = section->_parent_element_num; } } /* Construction of nodal connectivities */ /*---------------------------------------*/ for (cell_counter = 0 ; cell_counter < n_extr_cells ; cell_counter++) { if (extr_cells != NULL) cell_id = extr_cells[cell_counter] - 1; else cell_id = cell_counter; cell_type = _nodal_cell_from_desc(cell_id, n_face_lists, face_list_shift, face_vertex_idx, face_vertex_num, cell_face_idx, cell_face_num, cell_vtx_tria, cell_vtx_quad); section = sections[cell_type]; p_cell_vertex = section->_vertex_num + ( n_elements_type[cell_type] * fvm_nodal_n_vertices_element[cell_type]); switch (cell_type) { case FVM_CELL_TETRA: retcode = _nodal_from_desc_cnv_cel_tetra(cell_vtx_tria, p_cell_vertex); break; case FVM_CELL_PYRAM: retcode = _nodal_from_desc_cnv_cel_pyram(cell_vtx_tria, cell_vtx_quad, p_cell_vertex); break; case FVM_CELL_PRISM: retcode = _nodal_from_desc_cnv_cel_prism(cell_vtx_tria, cell_vtx_quad, p_cell_vertex); break; case FVM_CELL_HEXA: retcode = _nodal_from_desc_cnv_cel_hexa(cell_vtx_quad, p_cell_vertex); break; default: retcode = FVM_NODAL_FROM_DESC_SUCCESS; break; } /* Temporary value of parent cell num based on local cell id (so that the list of polyhedra given as an argument to _fvm_nodal_extract_polyhedra() below is based on the local numbering, like the cell->face connectivity */ section->_parent_element_num[n_elements_type[cell_type]] = cell_id + 1; n_elements_type[cell_type] += 1; if (retcode == FVM_NODAL_FROM_DESC_WORIENT) n_orient_pbs += 1; else if (retcode == FVM_NODAL_FROM_DESC_FAILURE) bft_error(__FILE__, __LINE__, 0, _("Incoherent connectivity for cell %d\n"), cell_id + 1); } fvm_parall_counter(&n_orient_pbs, 1); if (n_orient_pbs > 0) bft_printf(_("Warning: Possible nodal connectivity orientation\n" " problems for at least %d cells\n"), n_orient_pbs); /* Extraction of remaining polyhedra */ /*-----------------------------------*/ if (sections[FVM_CELL_POLY] != NULL) _fvm_nodal_extract_polyhedra (sections[FVM_CELL_POLY], n_elements_type[FVM_CELL_POLY], sections[FVM_CELL_POLY]->parent_element_num, n_face_lists, face_list_shift, face_vertex_idx, face_vertex_num, cell_face_idx, cell_face_num, cell_face_list); /* We can now base the final value of the parent cell number on the parent (and not local) numbering */ _raise_sections_parent_num(FVM_N_ELEMENT_TYPES, sections, parent_cell_num); /* Add sections to nodal mesh structure */ _optimize_sections_parent_num(FVM_N_ELEMENT_TYPES, sections); _fvm_nodal_add_sections(this_nodal, FVM_N_ELEMENT_TYPES, sections); } /*---------------------------------------------------------------------------- * Convert and add faces from an descending connectivity mesh to a nodal mesh. * * If the optional filter list extr_faces[] argument is non-NULL, faces * {extr_faces[0], extr_faces[1], extr_faces[n_extr_faces - 1]} are converted * and added to the nodal mesh. If this filter is set to NULL, faces * {1, 2, ..., n_extr_faces} are considered. * * In addition, an optional parent_face_num[] array may also be given, in * case the descending connectivity mesh definition is based on a temporary * subset of a parent mesh, (corresponding to the parent_face_num[] list, * using 1 to n numbering), and the final nodal mesh element parent numbering * should correspond to that parent mesh and not the temporary subset. * * parameters: * this_nodal <-> nodal mesh structure * n_extr_faces <-- count of faces to add * extr_faces <-- optional filter list of faces to extract (1 to n) * n_face_lists <-- number of face lists * face_list_shift <-- face list to common number index shifts; * size: n_face_lists * face_vertex_idx <-- face -> vertex indexes (per face list) * face_vertex_num <-- face -> vertex numbers (per face list) * parent_face_num <-- face -> parent face number (1 to n) if non-trivial * (i.e. if face definitions correspond to a subset * of the parent mesh), NULL otherwise. *----------------------------------------------------------------------------*/ void fvm_nodal_from_desc_add_faces(fvm_nodal_t *const this_nodal, const fvm_lnum_t n_extr_faces, const fvm_lnum_t extr_faces[], const int n_face_lists, const fvm_lnum_t face_list_shift[], const fvm_lnum_t *const face_vertex_idx[], const fvm_lnum_t *const face_vertex_num[], const fvm_lnum_t parent_face_num[]) { int type_id; fvm_lnum_t face_counter, face_id; fvm_element_t face_type; fvm_lnum_t *p_vertex_idx; fvm_lnum_t *p_vertex_num; fvm_lnum_t poly_connect_size; fvm_lnum_t n_face_vertices; fvm_lnum_t n_elements_type[FVM_N_ELEMENT_TYPES]; fvm_gnum_t n_g_elements_type[FVM_N_ELEMENT_TYPES]; fvm_nodal_section_t *section; fvm_nodal_section_t *sections[FVM_N_ELEMENT_TYPES]; /* Initialization */ for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) { n_elements_type[type_id] = 0; sections[type_id] = NULL; } poly_connect_size = 0; /* Compute connectivity type */ /*---------------------------*/ for (face_counter = 0 ; face_counter < n_extr_faces ; face_counter++) { if (extr_faces != NULL) face_id = extr_faces[face_counter] - 1; else face_id = face_counter; n_face_vertices = _nodal_face_from_desc_size(face_id, n_face_lists, face_list_shift, face_vertex_idx); switch (n_face_vertices) { case 3: face_type = FVM_FACE_TRIA; break; case 4: face_type = FVM_FACE_QUAD; break; default: face_type = FVM_FACE_POLY; poly_connect_size += n_face_vertices; break; } n_elements_type[face_type] += 1; } #if 0 && defined(DEBUG) && !defined(NDEBUG) bft_printf("dbg : fvm_nodal_faces_from_desc\n" "n_tria = %d\n" "n_quad = %d\n" "n_poly = %d\n", n_elements_type[FVM_FACE_TRIA], n_elements_type[FVM_FACE_QUAD], n_elements_type[FVM_FACE_POLY]); #endif /* Set dimensions (and reset local counters) */ for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) n_g_elements_type[type_id] = n_elements_type[type_id]; fvm_parall_counter(n_g_elements_type, FVM_N_ELEMENT_TYPES); for (face_type = FVM_FACE_TRIA ; face_type <= FVM_FACE_POLY ; face_type++) { if (n_g_elements_type[face_type] > 0) { sections[face_type] = fvm_nodal_section_create(face_type); sections[face_type]->n_elements = n_elements_type[face_type]; this_nodal->n_faces += n_elements_type[face_type]; } n_elements_type[face_type] = 0; } /* Main memory allocations */ for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) { section = sections[type_id]; if (section != NULL) { if (section->type != FVM_FACE_POLY) { section->stride = fvm_nodal_n_vertices_element[type_id]; section->connectivity_size = section->stride * section->n_elements; BFT_MALLOC(section->_vertex_num, section->connectivity_size, fvm_lnum_t); section->vertex_num = section->_vertex_num; } else { section->stride = fvm_nodal_n_vertices_element[type_id]; section->connectivity_size = poly_connect_size; BFT_MALLOC(section->_vertex_index, section->n_elements + 1, fvm_lnum_t); BFT_MALLOC(section->_vertex_num, section->connectivity_size, fvm_lnum_t); section->vertex_index = section->_vertex_index; section->vertex_num = section->_vertex_num; section->_vertex_index[0] = 0; } } } for (type_id = 0 ; type_id < FVM_N_ELEMENT_TYPES ; type_id++) { section = sections[type_id]; if (section != NULL) { BFT_MALLOC(section->_parent_element_num, section->n_elements, fvm_lnum_t); section->parent_element_num = section->_parent_element_num; } } /* Construction of nodal connectivities */ /*---------------------------------------*/ for (face_counter = 0 ; face_counter < n_extr_faces ; face_counter++) { if (extr_faces != NULL) face_id = extr_faces[face_counter] - 1; else face_id = face_counter; n_face_vertices = _nodal_face_from_desc_size(face_id, n_face_lists, face_list_shift, face_vertex_idx); switch (n_face_vertices) { case 3: face_type = FVM_FACE_TRIA; section = sections[face_type]; p_vertex_num = section->_vertex_num + (n_elements_type[face_type] * 3); break; case 4: face_type = FVM_FACE_QUAD; section = sections[face_type]; p_vertex_num = section->_vertex_num + (n_elements_type[face_type] * 4); break; default: face_type = FVM_FACE_POLY; section = sections[face_type]; p_vertex_idx = section->_vertex_index + n_elements_type[face_type]; *(p_vertex_idx + 1) = *p_vertex_idx + n_face_vertices; p_vertex_num = section->_vertex_num + (*p_vertex_idx); break; } _nodal_face_from_desc_copy(face_id, n_face_lists, face_list_shift, face_vertex_idx, face_vertex_num, p_vertex_num); section->_parent_element_num[n_elements_type[face_type]] = face_id + 1; n_elements_type[face_type] += 1; } /* We can now base the final value of the parent face number on the parent (and not local) numbering */ _raise_sections_parent_num(FVM_N_ELEMENT_TYPES, sections, parent_face_num); /* Add sections to nodal mesh structure */ _optimize_sections_parent_num(FVM_N_ELEMENT_TYPES, sections); _fvm_nodal_add_sections(this_nodal, FVM_N_ELEMENT_TYPES, sections); } /*----------------------------------------------------------------------------*/ #ifdef __cplusplus } #endif /* __cplusplus */