/*============================================================================ * Triangulation of nodal mesh sections *============================================================================*/ /* 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 #include /*---------------------------------------------------------------------------- * BFT library headers *----------------------------------------------------------------------------*/ #include #include /*---------------------------------------------------------------------------- * Local headers *----------------------------------------------------------------------------*/ #include #include #include #include #include #include /*---------------------------------------------------------------------------- * Header for the current file *----------------------------------------------------------------------------*/ #include /*----------------------------------------------------------------------------*/ #ifdef __cplusplus extern "C" { #if 0 } /* Fake brace to force Emacs auto-indentation back to column 0 */ #endif #endif /* __cplusplus */ /*============================================================================ * Static global variables *============================================================================*/ /*============================================================================ * Local macro definitions *============================================================================*/ /*============================================================================ * Private function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Create a new mesh section from projection of a given extruded section * to its base plane. * * The algorithm used here always leads to one edge per face. * The global numbering of the new section is not generated here, as it may * simply be transferred from the initial section. * * parameters: * dim <-- spatial dimension * vertex_coords <-- associated vertex coordinates array * parent_vertex_num <-- optional indirection to vertex coordinates * base_section <-- pointer to structure that should be triangulated * error_count --> number of triangulation errors counter (optional) * * returns: * pointer to created nodal mesh section representation structure *----------------------------------------------------------------------------*/ static fvm_nodal_section_t * _faces_to_edges(int dim, int chosen_axis, const fvm_coord_t vertex_coords[], const fvm_lnum_t parent_vertex_num[], const fvm_nodal_section_t *base_section, _Bool *selected_vertices, fvm_lnum_t *error_count) { fvm_lnum_t n_vertices, n_elements; fvm_lnum_t i, j, k, vertex_id; fvm_nodal_section_t *ret_section = NULL; /* Initialization */ if (error_count != NULL) *error_count = 0; n_elements = base_section->n_elements; /* Create new section */ /*--------------------*/ ret_section = fvm_nodal_section_create(FVM_EDGE); ret_section->n_elements = base_section->n_elements; ret_section->stride = 2; ret_section->connectivity_size = ret_section->stride * ret_section->n_elements; BFT_MALLOC(ret_section->_vertex_num, ret_section->connectivity_size, fvm_lnum_t); ret_section->vertex_num = ret_section->_vertex_num; if (base_section->parent_element_num != NULL) { BFT_MALLOC(ret_section->_parent_element_num, ret_section->n_elements, fvm_lnum_t); ret_section->parent_element_num = ret_section->_parent_element_num; } /* Main loop on section face elements */ /*------------------------------------*/ for (i = 0 ; i < n_elements ; i++) { fvm_lnum_t tmp_id[2], tmp_selected_edge[2]; fvm_coord_t tmp_edge_center; fvm_lnum_t selected_edge[2] = {0, 0}; fvm_coord_t selected_edge_center = 0; /* Compute number of vertices on each face */ if (base_section->vertex_index != NULL) { n_vertices = base_section->vertex_index[i+1] - base_section->vertex_index[i]; vertex_id = base_section->vertex_index[i]; } else { n_vertices = base_section->stride; vertex_id = base_section->stride * i; } /* Initialize the edge selection */ selected_edge[0] = base_section->vertex_num[vertex_id + n_vertices - 1]; selected_edge[1] = base_section->vertex_num[vertex_id]; if (parent_vertex_num == NULL) selected_edge_center = 0.5 * ( vertex_coords[(selected_edge[0]-1)*dim + chosen_axis] + vertex_coords[(selected_edge[1]-1)*dim + chosen_axis] ); else { tmp_id[0] = parent_vertex_num[selected_edge[0] - 1]; tmp_id[1] = parent_vertex_num[selected_edge[1] - 1]; selected_edge_center = 0.5 * ( vertex_coords[(tmp_id[0]-1) * dim + chosen_axis] + vertex_coords[(tmp_id[1]-1) * dim + chosen_axis] ); } /* Loop on face's vertices */ for (k = 1; k < n_vertices ; k++) { tmp_selected_edge[0] = base_section->vertex_num[vertex_id + k - 1]; tmp_selected_edge[1] = base_section->vertex_num[vertex_id + k]; if (parent_vertex_num == NULL) tmp_edge_center = 0.5 * ( vertex_coords[(tmp_selected_edge[0]-1)*dim + chosen_axis] + vertex_coords[(tmp_selected_edge[1]-1)*dim + chosen_axis]); else { tmp_id[0] = parent_vertex_num[tmp_selected_edge[0] - 1]; tmp_id[1] = parent_vertex_num[tmp_selected_edge[1] - 1]; tmp_edge_center = 0.5 * ( vertex_coords[(tmp_id[0]-1)*dim + chosen_axis] + vertex_coords[(tmp_id[1]-1)*dim + chosen_axis]); } if (tmp_edge_center < selected_edge_center) { selected_edge_center = tmp_edge_center; selected_edge[0] = tmp_selected_edge[0]; selected_edge[1] = tmp_selected_edge[1]; } } /* End of loop on element's vertices */ selected_vertices[selected_edge[0]-1] = true; selected_vertices[selected_edge[1]-1] = true; /* Define vertex_num */ for (j = 0; j < 2; j++) ret_section->_vertex_num[i*2 + j] = selected_edge[j]; if (base_section->parent_element_num != NULL) ret_section->_parent_element_num[i] = base_section->parent_element_num[i]; } /* End of loop on elements of the section */ /* Return new section */ return ret_section; } /*---------------------------------------------------------------------------- * Delete unused vertices after extracting edges from faces. * * parameters: * this_nodal <-> nodal mesh structure * selected_vertices <-- number of triangulation errors counter (optional) * * returns: * pointer to created nodal mesh section representation structure *----------------------------------------------------------------------------*/ static void _compact_mesh(fvm_nodal_t *this_nodal, _Bool *selected_vertices) { int i, j, i_section; int idx = 0; fvm_lnum_t *old_to_new = NULL; fvm_lnum_t *new_to_old = NULL; fvm_lnum_t *new_parent_vtx_num = NULL; fvm_coord_t *new_coords = NULL; fvm_lnum_t n_vertices_after = 0; fvm_lnum_t n_vertices_before = this_nodal->n_vertices; const int dim = this_nodal->dim; /* Count the new number of vertices */ for (i = 0; i < n_vertices_before; i++) { if (selected_vertices[i] == true) idx++; } n_vertices_after = idx; assert(n_vertices_after > 0); if (n_vertices_after == n_vertices_before) return; BFT_MALLOC(new_to_old, n_vertices_after, fvm_lnum_t); BFT_MALLOC(old_to_new, n_vertices_before, fvm_lnum_t); for (i = 0, idx = 0; i < n_vertices_before; i++) { old_to_new[i] = -1; if (selected_vertices[i] == true) { new_to_old[idx] = i + 1; old_to_new[i] = idx + 1; idx++; } } /* Adapt nodal structure if required */ /* --------------------------------- */ if (this_nodal->_vertex_coords != NULL) { /* generate a new coordinate array */ BFT_MALLOC(new_coords, dim * n_vertices_after, fvm_coord_t); if (this_nodal->_parent_vertex_num != NULL) { BFT_FREE(this_nodal->_parent_vertex_num); this_nodal->parent_vertex_num = NULL; } for (i = 0, idx = 0; i < n_vertices_before; i++) { if (selected_vertices[i] == true) { for (j = 0; j < dim; j++) new_coords[idx*dim + j] = this_nodal->vertex_coords[i*dim + j]; } } } /* End if _coords != NULL */ else { /* this_nodal->_coords == NULL */ if (this_nodal->parent_vertex_num != NULL) { /* generate a new parent_vertex_num */ BFT_MALLOC(new_parent_vtx_num, n_vertices_after, fvm_lnum_t); for (i = 0, idx = 0; i < n_vertices_before; i++) { if (selected_vertices[i] == true) new_parent_vtx_num[idx++] = this_nodal->parent_vertex_num[i]; } if (this_nodal->_parent_vertex_num != NULL) BFT_FREE(this_nodal->_parent_vertex_num); this_nodal->_parent_vertex_num = new_parent_vtx_num; this_nodal->parent_vertex_num = this_nodal->_parent_vertex_num; } /* End if parent_vertex_num != NULL */ } /* End if this_nodal->_coords == NULL */ /* Adapt _vertex_num in each section */ /* --------------------------------- */ for (i_section = 0; i_section < this_nodal->n_sections; i_section++) { fvm_nodal_section_t *section = this_nodal->sections[i_section]; const fvm_lnum_t n_connect = section->stride * section->n_elements; if (section->type == FVM_EDGE) { if (section->_vertex_num == NULL) BFT_MALLOC(section->_vertex_num, n_connect, fvm_lnum_t); for (j = 0; j < n_connect; j++) section->_vertex_num[j] = old_to_new[section->vertex_num[j] - 1]; section->vertex_num = section->_vertex_num; } } /* Adapt global_vertex_num if required */ /* ----------------------------------- */ if (this_nodal->global_vertex_num != NULL) { const fvm_gnum_t *old_global_vtx_num = fvm_io_num_get_global_num(this_nodal->global_vertex_num); fvm_io_num_t *new_global_vtx_num = fvm_io_num_create(new_to_old, old_global_vtx_num, n_vertices_after, 0); fvm_io_num_destroy(this_nodal->global_vertex_num); this_nodal->global_vertex_num = new_global_vtx_num; } this_nodal->n_vertices = n_vertices_after; } /*============================================================================ * Public function definitions *============================================================================*/ /*---------------------------------------------------------------------------- * Project an extruded mesh to its base plane. * * parameters: * this_nodal <-> pointer to structure that should be cut in edges * chosen_axis <-- indicate which axis is selected to extract edges * error_count --> number of triangulation errors counter (optional) *----------------------------------------------------------------------------*/ void fvm_nodal_project(fvm_nodal_t *this_nodal, int chosen_axis, fvm_lnum_t *error_count) { int i; fvm_lnum_t n_vertices = 0; fvm_lnum_t n_edges = 0; fvm_lnum_t section_error_count = 0; _Bool *selected_vertices = NULL; assert(this_nodal != NULL); n_vertices = this_nodal->n_vertices; BFT_MALLOC(selected_vertices, n_vertices, _Bool); for (i = 0; i < n_vertices; i++) selected_vertices[i] = false; /* Now extract edges and update new section list */ for (i = 0; i < this_nodal->n_sections; i++) { fvm_nodal_section_t *e_section = NULL; fvm_nodal_section_t *_section = this_nodal->sections[i]; if (_section->entity_dim == 2) { e_section = _faces_to_edges(this_nodal->dim, chosen_axis, this_nodal->vertex_coords, this_nodal->parent_vertex_num, _section, selected_vertices, §ion_error_count); if (error_count != NULL) *error_count += section_error_count; /* Transfer global numbering if present */ if (_section->global_element_num != NULL) { e_section->global_element_num = _section->global_element_num; _section->global_element_num = NULL; } fvm_nodal_section_destroy(_section); this_nodal->sections[i] = e_section; n_edges += e_section->n_elements; } } /* Adapt parent_vertex_num, coords and global_vertex_num */ _compact_mesh(this_nodal, selected_vertices); this_nodal->n_faces = 0; this_nodal->n_edges = n_edges; BFT_FREE(selected_vertices); } /*----------------------------------------------------------------------------*/ #ifdef __cplusplus } #endif /* __cplusplus */