/*============================================================================
* 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 <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/*----------------------------------------------------------------------------
* BFT library headers
*----------------------------------------------------------------------------*/
#include <bft_mem.h>
#include <bft_printf.h>
/*----------------------------------------------------------------------------
* Local headers
*----------------------------------------------------------------------------*/
#include <fvm_config_defs.h>
#include <fvm_defs.h>
#include <fvm_io_num.h>
#include <fvm_nodal.h>
#include <fvm_nodal_priv.h>
#include <fvm_nodal_project.h>
/*----------------------------------------------------------------------------
* Header for the current file
*----------------------------------------------------------------------------*/
#include <fvm_nodal_triangulate.h>
/*----------------------------------------------------------------------------*/
#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 */
syntax highlighted by Code2HTML, v. 0.9.1