/*============================================================================
* 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_triangulate.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
*============================================================================*/
#define _CROSS_PRODUCT_3D(cross_v1_v2, v1, v2) ( \
cross_v1_v2[0] = v1[1]*v2[2] - v1[2]*v2[1], \
cross_v1_v2[1] = v1[2]*v2[0] - v1[0]*v2[2], \
cross_v1_v2[2] = v1[0]*v2[1] - v1[1]*v2[0] )
#define _DOT_PRODUCT_3D(v1, v2) ( \
v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])
/*============================================================================
* Private function definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* Triangulate a quadrangle.
*
* A convex quadrangle is divided into two triangles along its shortest
* diagonal. A non-convex quadrangle may only be divided along the diagonal
* which lies inside the quadrangle.
*
* If the quadrangle_vertices argument is NULL, 1, 2, ...,n local numbering
* is implied.
*
* parameters:
* dim <-- spatial dimension (2 or 3).
* coords <-- coordinates of the triangulation's vertices.
* parent_vertex_num <-- optional indirection to vertex coordinates
* quadrangle_vertices <-- polygon connectivity; size: n_vertices or empty.
* triangle_vertices --> triangles connectivity; size: 2 * 3.
*
* returns:
* number of resulting triangles.
*----------------------------------------------------------------------------*/
static int
_triangulate_quadrangle(int dim,
const fvm_coord_t coords[],
const fvm_lnum_t parent_vertex_num[],
const fvm_lnum_t quadrangle_vertices[],
fvm_lnum_t triangle_vertices[])
{
int i, j;
double d2_02, d2_13;
int o_count = 0, o_id = 0;
fvm_lnum_t vertex_id[4] = {0, 1, 2, 3};
double v1[3] = {0.0, 0.0, 0.0}, v2[3] = {0.0, 0.0, 0.0};
double n0[3] = {0.0, 0.0, 0.0}, ni[3] = {0.0, 0.0, 0.0};
if (quadrangle_vertices != NULL) {
for (i = 0; i < 4 ; i++)
vertex_id[i] = quadrangle_vertices[i] - 1;
}
if (parent_vertex_num != NULL) {
for (i = 0; i < 4 ; i++)
vertex_id[i] = parent_vertex_num[i] - 1;
}
/* Check for an obtuse angle */
for (i = 0; i < dim; i++) {
v1[i] = coords[vertex_id[1]*dim + i] - coords[vertex_id[0]*dim + i];
v2[i] = coords[vertex_id[3]*dim + i] - coords[vertex_id[0]*dim + i];
}
_CROSS_PRODUCT_3D(n0, v1, v2);
for (j = 1; j < 4; j++) {
for (i = 0; i < dim; i++) {
v1[i] = coords[vertex_id[(j+1)%4]*dim + i] - coords[vertex_id[j]*dim + i];
v2[i] = coords[vertex_id[ j-1 ]*dim + i] - coords[vertex_id[0]*dim + i];
}
_CROSS_PRODUCT_3D(ni, v1, v2);
if (_DOT_PRODUCT_3D(n0, ni) < 0) {
o_count++;
o_id = j;
}
}
/* With an obtuse angle, only one diagonal lies inside the quadrangle;
we define it as "shorter" */
if (o_count > 0) {
if (o_count > 1)
o_id = 0;
if (o_id%2 == 0) {
d2_02 = 0.;
d2_13 = 1.;
}
else {
d2_02 = 1.;
d2_13 = 0.;
}
}
/* With no obtuse angle, we choose the true shortest diagonal */
else {
for (i = 0; i < dim; i++) {
v1[i] = coords[vertex_id[2]*dim + i] - coords[vertex_id[0]*dim + i];
v2[i] = coords[vertex_id[3]*dim + i] - coords[vertex_id[1]*dim + i];
}
/* Now compute diagonal lengths (squared) */
d2_02 = v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2];
d2_13 = v2[0]*v2[0] + v2[1]*v2[1] + v2[2]*v2[2];
}
/* Now define triangulation */
if (quadrangle_vertices != NULL) {
if (d2_02 < d2_13) {
triangle_vertices[0] = quadrangle_vertices[0]; /* 1st triangle */
triangle_vertices[1] = quadrangle_vertices[1];
triangle_vertices[2] = quadrangle_vertices[2];
triangle_vertices[3] = quadrangle_vertices[2]; /* 2nd triangle */
triangle_vertices[4] = quadrangle_vertices[3];
triangle_vertices[5] = quadrangle_vertices[0];
}
else {
triangle_vertices[0] = quadrangle_vertices[0]; /* 1st triangle */
triangle_vertices[1] = quadrangle_vertices[1];
triangle_vertices[2] = quadrangle_vertices[3];
triangle_vertices[3] = quadrangle_vertices[2]; /* 2nd triangle */
triangle_vertices[4] = quadrangle_vertices[3];
triangle_vertices[5] = quadrangle_vertices[1];
}
}
else { /* if (quadrangle_vertices == NULL) */
if (d2_02 < d2_13) {
triangle_vertices[0] = 1; /* 1st triangle */
triangle_vertices[1] = 2;
triangle_vertices[2] = 3;
triangle_vertices[3] = 3; /* 2nd triangle */
triangle_vertices[4] = 4;
triangle_vertices[5] = 1;
}
else {
triangle_vertices[0] = 1; /* 1st triangle */
triangle_vertices[1] = 2;
triangle_vertices[2] = 4;
triangle_vertices[3] = 3; /* 2nd triangle */
triangle_vertices[4] = 4;
triangle_vertices[5] = 2;
}
}
/* Return number of triangles (for consistency with polygon triangulation) */
return 2;
}
/*----------------------------------------------------------------------------
* Create a new mesh section from the triangulation of a given 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 *
_triangulate_section(int dim,
const fvm_coord_t vertex_coords[],
const fvm_lnum_t parent_vertex_num[],
const fvm_nodal_section_t *base_section,
fvm_lnum_t *error_count)
{
fvm_lnum_t n_vertices, n_triangles, n_elements;
fvm_lnum_t n_vertices_max, n_triangles_tot;
fvm_lnum_t i, j, k, triangle_id, vertex_id;
fvm_triangulate_state_t *state = NULL;
fvm_lnum_t *n_sub_elements = NULL;
fvm_nodal_section_t *ret_section = NULL;
/* Initialization */
if (error_count != NULL)
*error_count = 0;
n_elements = base_section->n_elements;
if (base_section->global_element_num != NULL)
BFT_MALLOC(n_sub_elements, n_elements, fvm_lnum_t);
/* Count expected total and local numbers of triangles */
n_vertices_max = 0;
n_triangles_tot = 0;
if (base_section->vertex_index != NULL) {
for (i = 0 ; i < n_elements ; i++) {
n_vertices = base_section->vertex_index[i+1]
- base_section->vertex_index[i];
n_triangles_tot += n_vertices - 2;
if (n_vertices > n_vertices_max)
n_vertices_max = n_vertices;
}
}
else if (base_section->stride == 4) {
n_triangles_tot = base_section->n_elements * 2;
n_vertices_max = 4;
}
else if (base_section->stride == 3) {
n_triangles_tot = base_section->n_elements;
n_vertices_max = 3;
}
/* Allocate memory and state variables */
if (n_vertices_max > 4 && base_section->vertex_index != NULL)
state = fvm_triangulate_state_create(n_vertices_max);
/* Create new section */
/*--------------------*/
ret_section = fvm_nodal_section_create(FVM_FACE_TRIA);
ret_section->n_elements = n_triangles_tot;
ret_section->stride = 3;
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;
BFT_MALLOC(ret_section->_parent_element_num,
ret_section->n_elements,
fvm_lnum_t);
ret_section->parent_element_num = ret_section->_parent_element_num;
triangle_id = 0;
/* Main loop on section face elements */
/*------------------------------------*/
for (i = 0 ; i < n_elements ; i++) {
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;
}
n_triangles = 0;
/* If face must be subdivided */
if (n_vertices >= 4) {
if (n_vertices == 4)
n_triangles = _triangulate_quadrangle(dim,
vertex_coords,
parent_vertex_num,
( base_section->vertex_num
+ vertex_id),
( ret_section->_vertex_num
+ triangle_id*3));
else {
n_triangles = fvm_triangulate_polygon(dim,
n_vertices,
vertex_coords,
parent_vertex_num,
( base_section->vertex_num
+ vertex_id),
( ret_section->_vertex_num
+ triangle_id*3),
state);
if (n_triangles != (n_vertices - 2) && error_count != NULL)
*error_count += 1;
}
if (base_section->parent_element_num != NULL) {
for (j = 0; j < n_triangles; j++)
ret_section->_parent_element_num[triangle_id + j]
= base_section->parent_element_num[i];
}
else {
for (j = 0; j < n_triangles; j++)
ret_section->_parent_element_num[triangle_id + j] = i+1;
}
triangle_id += n_triangles;
}
/* Otherwise, face must simply be copied */
else if (n_vertices == 3) {
n_triangles = 1;
for (k = 0; k < 3; k++)
ret_section->_vertex_num[triangle_id*3 + k]
= base_section->vertex_num[i*3 + k];
if (base_section->parent_element_num != NULL)
ret_section->_parent_element_num[triangle_id]
= base_section->parent_element_num[i];
else
ret_section->_parent_element_num[triangle_id] = i+1;
triangle_id += 1;
}
if (n_sub_elements != NULL)
n_sub_elements[i] = n_triangles;
}
/* Free memory and state variables */
if (n_vertices_max > 4 && base_section->vertex_index != NULL)
state = fvm_triangulate_state_destroy(state);
/* Now update global numbering if present */
if (base_section->global_element_num != NULL) {
ret_section->global_element_num
= fvm_io_num_create_from_sub(base_section->global_element_num,
n_sub_elements);
if (n_sub_elements != NULL)
BFT_FREE(n_sub_elements);
}
/* Return new (triangulated) section */
return ret_section;
}
/*----------------------------------------------------------------------------
* Create new mesh sections from the triangulation of a given 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
* new_sections --> array of triangle and quadrangle sections
* error_count --> number of triangulation errors counter (optional)
*----------------------------------------------------------------------------*/
static void
_triangulate_section_polygons(int dim,
const fvm_coord_t vertex_coords[],
const fvm_lnum_t parent_vertex_num[],
const fvm_nodal_section_t *base_section,
fvm_nodal_section_t *new_sections[2],
fvm_lnum_t *error_count)
{
int type_id;
fvm_lnum_t n_vertices, n_triangles, n_elements;
fvm_lnum_t n_vertices_max, n_triangles_max;
fvm_lnum_t i, j, k, vertex_id;
int *triangle_vertices = NULL;
fvm_element_t element_type[2] = {FVM_FACE_TRIA,
FVM_FACE_QUAD};
fvm_lnum_t n_elements_tot[2] = {0, 0}; /* New triangles/quadrangles */
fvm_lnum_t element_id[2] = {0, 0};
fvm_lnum_t *n_sub_elements[2] = {NULL, NULL};
fvm_triangulate_state_t *state = NULL;
/* Initialization */
if (error_count != NULL)
*error_count = 0;
new_sections[0] = NULL, new_sections[1] = NULL;
n_elements = base_section->n_elements;
/* Count expected total and local numbers of triangles */
n_vertices_max = 0;
n_triangles_max = 0;
if (base_section->vertex_index != NULL) {
for (i = 0 ; i < n_elements ; i++) {
n_vertices = base_section->vertex_index[i+1]
- base_section->vertex_index[i];
if (n_vertices == 4)
n_elements_tot[1] += 1;
else
n_elements_tot[0] += n_vertices - 2;
if (n_vertices > n_vertices_max)
n_vertices_max = n_vertices;
}
}
else
return;
n_triangles_max = n_vertices_max - 2;
/* Allocate memory and state variables */
if (n_vertices_max > 4 && base_section->vertex_index != NULL) {
BFT_MALLOC(triangle_vertices, n_triangles_max*3, int);
state = fvm_triangulate_state_create(n_vertices_max);
}
/* Create new sections */
/*---------------------*/
for (type_id = 0; type_id < 2; type_id++) {
if (n_elements_tot[type_id] > 0) {
fvm_nodal_section_t *_section;
_section = fvm_nodal_section_create(element_type[type_id]);
_section->n_elements = n_elements_tot[type_id];
_section->stride = fvm_nodal_n_vertices_element[element_type[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;
BFT_MALLOC(_section->_parent_element_num, _section->n_elements, fvm_lnum_t);
_section->parent_element_num = _section->_parent_element_num;
new_sections[type_id] = _section;
if (base_section->global_element_num != NULL)
BFT_MALLOC(n_sub_elements[type_id], n_elements, fvm_lnum_t);
}
}
/* Main loop on section face elements */
/*------------------------------------*/
for (i = 0 ; i < n_elements ; i++) {
fvm_nodal_section_t *_section;
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;
}
if (n_vertices == 4) {
type_id = 1;
}
else {
type_id = 0;
}
_section = new_sections[type_id];
/* If face must be subdivided */
if (n_vertices > 4) {
n_triangles = fvm_triangulate_polygon(dim,
n_vertices,
vertex_coords,
parent_vertex_num,
( base_section->vertex_num
+ vertex_id),
( _section->_vertex_num
+ element_id[0]*3),
state);
if (n_triangles != (n_vertices - 2) && error_count != NULL)
*error_count += 1;
if (base_section->parent_element_num != NULL) {
for (j = 0; j < n_triangles; j++)
_section->_parent_element_num[element_id[0] + j]
= base_section->parent_element_num[i];
}
else {
for (j = 0; j < n_triangles; j++)
_section->_parent_element_num[element_id[0] + j] = i+1;
}
element_id[0] += n_triangles;
/* Update sub-element counts */
if (n_sub_elements[0] != NULL)
(n_sub_elements[0])[i] = n_triangles;
if (n_sub_elements[1] != NULL)
(n_sub_elements[1])[i] = 0;
}
/* Otherwise, face must simply be copied */
else {
for (k = 0; k < _section->stride; k++)
_section->_vertex_num[element_id[type_id]*_section->stride + k]
= base_section->vertex_num[i*_section->stride + k];
if (base_section->parent_element_num != NULL)
_section->_parent_element_num[element_id[type_id]]
= base_section->parent_element_num[i];
else {
_section->_parent_element_num[element_id[type_id]] = i+1;
}
element_id[type_id] += 1;
/* Update sub-element counts */
if (n_sub_elements[type_id] != NULL)
(n_sub_elements[type_id])[i] = 1;
if (n_sub_elements[(type_id+1)%2] != NULL)
(n_sub_elements[(type_id+1)%2])[i] = 0;
}
}
/* Free memory and state variables */
if (n_vertices_max > 4 && base_section->vertex_index != NULL) {
BFT_FREE(triangle_vertices);
state = fvm_triangulate_state_destroy(state);
}
/* Now update global numbering if present */
for (type_id = 0; type_id < 2; type_id++) {
if (n_sub_elements[type_id] != NULL) {
(new_sections[type_id])->global_element_num
= fvm_io_num_create_from_sub(base_section->global_element_num,
n_sub_elements[type_id]);
}
BFT_FREE(n_sub_elements[type_id]);
}
}
/*============================================================================
* Public function definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* Triangulate all sections of a nodal mesh.
*
* parameters:
* this_nodal <-> pointer to structure that should be triangulated
* error_count --> number of triangulation errors counter (optional)
*----------------------------------------------------------------------------*/
void
fvm_nodal_triangulate(fvm_nodal_t *this_nodal,
fvm_lnum_t *error_count)
{
int i;
fvm_lnum_t section_error_count = 0;
fvm_lnum_t n_faces = 0;
assert(this_nodal != NULL);
/* Now triangulate and update new section list */
for (i = 0; i < this_nodal->n_sections; i++) {
fvm_nodal_section_t *t_section;
fvm_nodal_section_t *_section = this_nodal->sections[i];
if (_section->entity_dim == 2) {
if (_section->type != FVM_FACE_TRIA) {
t_section = _triangulate_section(this_nodal->dim,
this_nodal->vertex_coords,
this_nodal->parent_vertex_num,
_section,
§ion_error_count);
if (error_count != NULL)
*error_count += section_error_count;
fvm_nodal_section_destroy(_section);
this_nodal->sections[i] = t_section;
n_faces += t_section->n_elements;
}
else
n_faces +=_section->n_elements;
}
}
this_nodal->n_faces = n_faces;
}
/*----------------------------------------------------------------------------
* Triangulate polygonal sections of a nodal mesh.
*
* parameters:
* this_nodal <-> pointer to structure that should be triangulated
* error_count --> number of triangulation errors counter (optional)
*----------------------------------------------------------------------------*/
void
fvm_nodal_triangulate_polygons(fvm_nodal_t *this_nodal,
fvm_lnum_t *error_count)
{
int i, j;
fvm_nodal_section_t *new_sections[2];
int n_sections = 0;
fvm_lnum_t section_error_count = 0;
fvm_lnum_t n_faces = 0;
fvm_nodal_section_t **sections = NULL;
assert(this_nodal != NULL);
/* Best estimate for new section list size: polygonal sections
may lead to 2 sections (triangles + quads) */
for (i = 0; i < this_nodal->n_sections; i++) {
fvm_nodal_section_t *_section = this_nodal->sections[i];
if (_section->type == FVM_FACE_POLY)
n_sections += 2;
else
n_sections += 1;
}
BFT_MALLOC(sections, n_sections, fvm_nodal_section_t *);
/* Now triangulate and update new section list */
n_sections = 0;
for (i = 0; i < this_nodal->n_sections; i++) {
fvm_nodal_section_t *_section = this_nodal->sections[i];
if (_section->type == FVM_FACE_POLY) {
_triangulate_section_polygons(this_nodal->dim,
this_nodal->vertex_coords,
this_nodal->parent_vertex_num,
_section,
new_sections,
§ion_error_count);
if (error_count != NULL)
*error_count += section_error_count;
fvm_nodal_section_destroy(_section);
for (j = 0; j < 2; j++) {
if (new_sections[j] != NULL) {
sections[n_sections] = new_sections[j];
n_faces += (new_sections[j])->n_elements;
n_sections += 1;
}
}
}
else {
if (_section->entity_dim == 2)
n_faces += _section->n_elements;
sections[n_sections] = _section;
n_sections += 1;
}
}
/* Now replace section list */
BFT_FREE(this_nodal->sections);
BFT_REALLOC(sections, n_sections, fvm_nodal_section_t *);
this_nodal->n_sections = n_sections;
this_nodal->sections = sections;
this_nodal->n_faces = n_faces;
}
/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
}
#endif /* __cplusplus */
syntax highlighted by Code2HTML, v. 0.9.1