/*============================================================================
* Triangulation of a polygon
*============================================================================*/
/*
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) 2005-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>
/*----------------------------------------------------------------------------
* Local headers
*----------------------------------------------------------------------------*/
#include <fvm_config_defs.h>
#include <fvm_defs.h>
/*----------------------------------------------------------------------------
* Header for the current file
*----------------------------------------------------------------------------*/
#include <fvm_triangulate.h>
/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
extern "C" {
#if 0
} /* Fake brace to force Emacs auto-indentation back to column 0 */
#endif
#endif /* __cplusplus */
/*============================================================================
* Local Type Definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* State of current triangulation (working structure)
*----------------------------------------------------------------------------*/
struct _fvm_triangulate_state_t {
int *triangle_vertices; /* current triangle vertices list */
fvm_coord_t *coords; /* vertex coordinates */
int *list_previous; /* indices of previous vertices in polygon;
size:n_vertices; */
int *list_next; /* indices of next vertices in polygon;
size:n_vertices; */
int *edge_vertices; /* edges connectivity; size: n_edges * 2 */
int *edge_neighbors; /* triangles sharing a given edge
(2 values per edge);
- (-1,-1): edge does not exist
- (x ,-1): boundary edge, triangle x
- (x ,-1): internal edge, triangles x
and y */
_Bool *edge_is_delaunay; /* Delaunay edge indicator */
_Bool *concave; /* Concave vertex indicator */
int n_vertices_max; /* Maximum number vertices; a larger
size requires resizing work arrays */
};
/*============================================================================
* 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])
#define _MODULE_3D(v) \
sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])
/*============================================================================
* Private function definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* Project a face in 3d on a plane parallel to this face, which is then
* treated as the (Oxy) plane.
*
* parameters:
* n_vertices <-- number of vertices defining the polygon.
* coords <-> coordinates of the polygon's vertices (3d in, 2d out)
*----------------------------------------------------------------------------*/
static void
_polygon_plane_3d(const int n_vertices,
fvm_coord_t coords[])
{
int i, j;
fvm_coord_t face_center[3], face_normal[3];
fvm_coord_t v1[3], v2[3], cross_v1_v2[3];
fvm_coord_t dot_v1_v2, module_v2;
fvm_coord_t tmp_coord;
double cost;
double sint;
#define _N_VERTICES_AUTO_MAX 20 /* Size of local temporary coordinates
buffer; above this size, allocation
is necessary */
fvm_coord_t _tmp_coords[_N_VERTICES_AUTO_MAX * 3];
fvm_coord_t *_tmp_coords_p = NULL;
fvm_coord_t *tmp_coords = _tmp_coords;
/*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/
/* Estimate position of polygon center */
for (i = 0; i < 3; i++)
face_center[i] = 0.;
for (i = 0; i < 3; i++) {
for (j = 0; j < n_vertices; j++)
face_center[i] += coords[j*3 + i];
face_center[i] /= n_vertices;
}
/* Estimate face normal */
for (i = 0; i < 3; i++)
face_normal[i] = 0.;
for (j = 0; j < n_vertices; j++) {
for (i = 0; i < 3; i++) {
v1[i] = coords[j*3 + i] - face_center[i];
if (j < n_vertices - 1)
v2[i] = coords[(j+1)*3 + i] - face_center[i];
else
v2[i] = coords[ i] - face_center[i];
}
face_normal[0] += v1[1]*v2[2] - v1[2]*v2[1];
face_normal[1] += v1[2]*v2[0] - v1[0]*v2[2];
face_normal[2] += v1[0]*v2[1] - v1[1]*v2[0];
}
/* Project coordinates in a plane parallel to face */
/*-------------------------------------------------*/
/* We place the coordinate system origin at the estimated face center */
for (j = 0; j < n_vertices; j++)
for (i = 0; i < 3; i++)
coords[j*3 + i] -= face_center[i];
if (FVM_ABS(face_normal[0]) > 1.e-12 || FVM_ABS(face_normal[1]) > 1.e-12) {
/* First rotation of axis (Oz) and angle (Ox, normal proj. on Oxy) */
if (n_vertices > _N_VERTICES_AUTO_MAX) {
BFT_MALLOC(_tmp_coords_p, n_vertices*3, fvm_coord_t);
tmp_coords = _tmp_coords_p;
}
v1[0] = 1.;
v1[1] = 0.;
v1[2] = 0.;
v2[0] = face_normal[0];
v2[1] = face_normal[1];
v2[2] = 0.;
_CROSS_PRODUCT_3D(cross_v1_v2, v1, v2);
dot_v1_v2 = _DOT_PRODUCT_3D(v1, v2);
module_v2 = _MODULE_3D(v2);
cost = dot_v1_v2 / module_v2;
if (cross_v1_v2[2] > 0.)
sint = _MODULE_3D(cross_v1_v2) / module_v2;
else
sint = - _MODULE_3D(cross_v1_v2) / module_v2;
for (j = 0; j < n_vertices; j++) {
tmp_coords[j*3 ] = cost*coords[j*3] + sint*coords[j*3 + 1];
tmp_coords[j*3 + 1] = -sint*coords[j*3] + cost*coords[j*3 + 1];
tmp_coords[j*3 + 2] = coords[j*3 +2];
}
/* Second rotation with axis (Oy) and angle (Oz', normal proj. on Ox'z) */
v1[0] = 0.;
v1[1] = 0.;
v1[2] = 1.;
v2[0] = sqrt(face_normal[0]*face_normal[0] + face_normal[1]*face_normal[1]);
v2[1] = 0.;
v2[2] = face_normal[2];
_CROSS_PRODUCT_3D(cross_v1_v2, v1, v2);
cost = _DOT_PRODUCT_3D(v1, v2) / _MODULE_3D(v2);
if (cross_v1_v2[2] > 0.)
sint = _MODULE_3D(cross_v1_v2) / _MODULE_3D(v2);
else
sint = -_MODULE_3D(cross_v1_v2) / _MODULE_3D(v2);
for (j = 0; j < n_vertices; j++) {
coords[j*3 ] = cost*tmp_coords[j*3] + sint*tmp_coords[j*3 + 2];
coords[j*3 + 1] = tmp_coords[j*3 + 1];
coords[j*3 + 2] = 0.;
}
if (_tmp_coords_p != NULL) {
BFT_FREE(_tmp_coords_p);
tmp_coords = NULL;
}
}
else {
/* We only need to set the vertices z coordinate to 0, possibly
swapping the coordinates in the (Oxy) projection plane. */
if (face_normal[2] > 0.)
for (j = 0; j < n_vertices; j++)
coords[j*3 + 2] = 0.;
else
for (j = 0; j < n_vertices; j++) {
tmp_coord = coords[j*3];
coords[j*3 ] = coords[j*3 + 1];
coords[j*3 + 1] = tmp_coord;
coords[j*3 + 2] = 0.;
}
}
/* Convert coords to 2d */
for (j = 0; j < n_vertices; j++) {
coords[j*2 ] = coords[j*3 ];
coords[j*2 + 1] = coords[j*3 + 1];
}
#undef _N_VERTICES_AUTO_MAX
}
/*----------------------------------------------------------------------------
* Check if a (2D projection of a) polygon vertex is convex.
*
* parameters:
* previous <-- index of previous vertex in polygon.
* current <-- index of current vertex in polygon.
* next <-- index of following vertex in polygon.
* coords <-- coordinates of the polygon's vertices (2d).
*----------------------------------------------------------------------------*/
static _Bool
_polygon_vertex_is_convex(const int previous,
const int current,
const int next,
const fvm_coord_t coords[])
{
/* sin(theta) = (v1 x v2) / (|v1| x |v2|), so the sign of the sine
is that of the cross product's z component (as v1 and v2 lie in
the same plane) */
if ( ( (coords[current*2] - coords[previous*2])
* (coords[next*2 + 1] - coords[current*2 + 1]))
- ( (coords[next*2] - coords[current*2])
* (coords[current*2 + 1] - coords[previous*2 + 1])) > 0.0)
return true;
else
return false;
}
/*----------------------------------------------------------------------------
* Check if a (2D projection of a) polygon vertex is an ear.
*
* parameters:
* n_vertices <-- number of vertices defining the polygon.
* current <-- index of current vertex in polygon.
* list_previous <-- indices of previous vertices in polygon.
* list_next <-- index of previous vertex in polygon.
* concave <-- flag concave vertices in polygon.
* coords <-- coordinates of the polygon's vertices (2d).
* epsilon <-- associated relative tolerance.
*----------------------------------------------------------------------------*/
static _Bool
_polygon_vertex_is_ear(const int n_vertices,
const int current,
const int list_previous[],
const int list_next[],
const _Bool concave[],
const fvm_coord_t coords[],
const double epsilon)
{
int i, previous, next;
double surf_2, x_iso, y_iso;
double vect1[2], vect2[2], vect3[2];
/* If no vertex is concave, we have an ear */
for (i = 0; i < n_vertices && concave[i] == false; i++);
if (i == n_vertices)
return true;
/* If current vertex is convex */
else {
if (concave[current] == false) {
/* We check if the triangle formed by
list_previous[current], current, and list_next[current]
contains a concave vertex */
previous = list_previous[current];
next = list_next[current];
vect2[0] = coords[current*2 ] - coords[previous*2 ];
vect2[1] = coords[current*2 + 1] - coords[previous*2 + 1];
vect3[0] = coords[next*2 ] - coords[previous*2 ];
vect3[1] = coords[next*2 + 1] - coords[previous*2 + 1];
surf_2 = vect2[0]*vect3[1] - vect3[0]*vect2[1];
for (i = list_next[next]; i != previous; i = list_next[i]) {
if (concave[i] == true) {
vect1[0] = coords[i*2 ] - coords[previous*2 ];
vect1[1] = coords[i*2 + 1] - coords[previous*2 + 1];
x_iso = (vect1[0]*vect3[1] - vect1[1]*vect3[0]) / surf_2;
y_iso = (vect2[0]*vect1[1] - vect2[1]*vect1[0]) / surf_2;
if ( (1.0 - x_iso - y_iso > - epsilon)
&& ( x_iso > - epsilon)
&& ( y_iso > - epsilon))
return false;
}
}
return true;
}
else
return false;
}
}
/*----------------------------------------------------------------------------
* Check if an edge (between two triangles) is locally Delaunay.
*
* We compute the power of a point compared to the circle circumscribed
* to the triangle formed by the three others (of which two define
* the diagonal considered).
*
* parameters:
* edge_vertex_0 <-- index of first edge vertex in coords.
* edge_vertex_1 <-- index of second edge vertex in coords.
* flip_vertex_0 <-- index of first flip edge vertex in coords.
* flip_vertex_1 <-- index of second flip edge vertex in coords.
* coords <-- coordinates of the triangulation's vertices (2d).
*----------------------------------------------------------------------------*/
static _Bool
_edge_is_locally_delaunay(const int edge_vertex_0,
const int edge_vertex_1,
const int flip_vertex_0,
const int flip_vertex_1,
const fvm_coord_t coords[])
{
double a, b, delta;
double lambda[4];
double x_center, y_center, radius;
double point_power;
lambda[0] = 2*(coords[edge_vertex_1*2 ] - coords[edge_vertex_0*2 ]);
lambda[1] = 2*(coords[edge_vertex_1*2 + 1] - coords[edge_vertex_0*2 + 1]);
lambda[2] = 2*(coords[flip_vertex_0*2 ] - coords[edge_vertex_0*2 ]);
lambda[3] = 2*(coords[flip_vertex_0*2 + 1] - coords[edge_vertex_0*2 + 1]);
delta = lambda[1]*lambda[2] - lambda[0]*lambda[3];
/* If the triangle is flat, we automatically switch diagonals to
avoid a division by zero. */
if (FVM_ABS(delta) < 1.e-12)
return false;
a = coords[edge_vertex_1*2 ] * coords[edge_vertex_1*2 ]
- coords[edge_vertex_0*2 ] * coords[edge_vertex_0*2 ]
+ coords[edge_vertex_1*2 + 1] * coords[edge_vertex_1*2 + 1]
- coords[edge_vertex_0*2 + 1] * coords[edge_vertex_0*2 + 1];
b = coords[flip_vertex_0*2 ] * coords[flip_vertex_0*2 ]
- coords[edge_vertex_0*2 ] * coords[edge_vertex_0*2 ]
+ coords[flip_vertex_0*2 + 1] * coords[flip_vertex_0*2 + 1]
- coords[edge_vertex_0*2 + 1] * coords[edge_vertex_0*2 + 1];
/* Center and radius of the circle passing through the vertices of
the diagonal and through a third vertex. */
x_center = (lambda[1]*b - lambda[3]*a)/delta;
y_center = (lambda[2]*a - lambda[0]*b)/delta;
radius = sqrt( ( (x_center - coords[edge_vertex_0*2 ])
* (x_center - coords[edge_vertex_0*2 ]))
+ ( (y_center - coords[edge_vertex_0*2 + 1])
* (y_center - coords[edge_vertex_0*2 + 1])));
/* Compute the power of the quadrilateral's fourth vertex compared
to the circle. */
point_power = ( (coords[flip_vertex_1*2 ] - x_center)
* (coords[flip_vertex_1*2 ] - x_center))
+ ( (coords[flip_vertex_1*2 + 1] - y_center)
* (coords[flip_vertex_1*2 + 1] - y_center))
- radius*radius;
/* Keep a slight margin in case there is no perceptible gain
in switching diagonals */
if (point_power > -1.e-12)
return true;
else
return false;
}
/*----------------------------------------------------------------------------
* Sort vertex indexes defining a triangle.
*
* parameters:
* triangle_vertices <-> triangles connectivity.
*----------------------------------------------------------------------------*/
static void
_triangle_by_sorted_vertices(int triangle_vertices[])
{
int i_min, i_max, i_mid;
/* First step */
if (triangle_vertices[0] < triangle_vertices[1]) {
i_min = triangle_vertices[0];
i_mid = triangle_vertices[1];
i_max = i_mid;
}
else {
i_min = triangle_vertices[1];
i_mid = triangle_vertices[0];
i_max = i_mid;
}
/* Second step */
if (triangle_vertices[2] < i_min) {
i_mid = i_min;
i_min = triangle_vertices[2];
}
else if (triangle_vertices[2] > i_max) {
i_max = triangle_vertices[2];
}
else {
i_mid = triangle_vertices[2];
}
/* Reordering */
triangle_vertices[0] = i_min;
triangle_vertices[1] = i_mid;
triangle_vertices[2] = i_max;
}
/*----------------------------------------------------------------------------
* Convert a given triangulation to a Delaunay triangulation using
* edge flips.
*
* parameters:
* n_vertices <-- number of vertices defining the triangulation.
* triangle_vertices <-> triangles connectivity.
* edge_vertices --> edges connectivity.
* edge_neigbors --> triangles sharing a given edge (2 values per edge);
* - (-1,-1): edge does not exist
* - (x ,-1): boundary edge, triangle x
* - (x ,-1): internal edge, triangles x and y
* edge_is_delaunay --> delaunay edge indicator.
* coords <-- coordinates of the triangulation's vertices (2d).
*----------------------------------------------------------------------------*/
static void
_polygon_delaunay_flip(const int n_vertices,
int triangle_vertices[],
int edge_vertices[],
int edge_neighbors[],
_Bool edge_is_delaunay[],
const fvm_coord_t coords[])
{
int triangle_id, edge_id, vertex_id;
int triangle_id_0, triangle_id_1;
int current_edge_id, flip_edge_id;
int vertex_flip[2];
int i, i_0, i_1, i_min, i_max, j;
_Bool face_is_delaunay, edge_locally_delaunay;
_Bool restart, convex_quad;
const int n_edges = n_vertices*(n_vertices - 1)/2;
const int n_triangles = n_vertices - 2;
int i_previous[2] = {-1, -1}, i_next[2] = {-1, -1};
/*
There are (n_vertices*(n_vertices - 1)/2) possible edge combinations.
An edge's number can be given by:
-> SUM(k = 0, k < i_min) {n_vertices - k} + (i_max - i_min)
i.e. n_vertices*i_min - i_min*(i_min+1)/2 + (i_max - i_min)
An edge's index equals it's number - 1.
Be careful with the order of arguments! arg[0]=i_min, arg[1]=i_max
*/
#undef _EDGE_INDEX
#define _EDGE_INDEX(i_min, i_max) \
(n_vertices*i_min - i_min*(i_min+1)/2 + i_max-i_min - 1)
/* Initialization */
/*----------------*/
for (i_0 = 0; i_0 < n_vertices; i_0++) {
for (i_1 = i_0 + 1; i_1 < n_vertices; i_1++) {
edge_id = _EDGE_INDEX(i_0, i_1);
/* Define edges */
edge_vertices[2*edge_id ] = i_0;
edge_vertices[2*edge_id + 1] = i_1;
/*
Liste of triangles sharing an edge:
- (-1,-1): edge does not exist
- (x ,-1): boundary edge, triangle x
- (x ,-1): internal edge, triangles x and y
*/
edge_neighbors[2*edge_id ] = -1;
edge_neighbors[2*edge_id + 1] = -1;
/* Initialize an array indicating if an edge is locally Delaunay */
edge_is_delaunay[edge_id] = true;
}
}
/* First traversal of triangles to determine initial neighborhood,
as well as the list of edges which are not locally Delaunay */
for (j = 0; j < n_triangles; j++) {
for (i = 0; i < 3; i++) {
i_0 = triangle_vertices[(j*3) + i ];
i_1 = triangle_vertices[(j*3) + ((i+1)%3)];
i_min = FVM_MIN(i_0, i_1);
i_max = FVM_MAX(i_0, i_1);
edge_id = _EDGE_INDEX(i_min, i_max);
/* Update edge neighbors */
if (edge_neighbors[2*edge_id] == -1)
edge_neighbors[2*edge_id] = j;
else
edge_neighbors[2*edge_id + 1] = j;
/* If edge is not on the boundary: */
if ( !(i_max == i_min + 1)
&& !(i_min == 0 && i_max == n_vertices - 1))
edge_is_delaunay[edge_id] = false;
}
}
/* Main flip algorithm */
/*---------------------*/
edge_id = 0;
restart = false;
face_is_delaunay = false;
while(face_is_delaunay == false) {
if (edge_is_delaunay[edge_id] == false) {
edge_is_delaunay[edge_id] = true;
i_0 = edge_vertices[2*edge_id];
i_1 = edge_vertices[2*edge_id + 1];
for (j = 0; j < 2; j++) { /* Loop on triangles on each side of edge */
triangle_id = edge_neighbors[2*edge_id + j];
for (i = 0; i < 3; i++) { /* Loop on triangle's vertices */
vertex_id = triangle_vertices[3*triangle_id + i];
/* Seek opposite vertices */
if (vertex_id != i_0 && vertex_id != i_1)
vertex_flip[j] = vertex_id;
/* Seek preceding and following vertices */
if ( vertex_id == i_0
&& triangle_vertices[3*triangle_id + ((i + 1)%3)] == i_1) {
i_previous[0] = triangle_vertices[3*triangle_id + ((i + 2)%3)];
i_next[1] = i_previous[0];
}
else if ( vertex_id == i_1
&& triangle_vertices[3*triangle_id + ((i + 1)%3)] == i_0) {
i_next[0] = triangle_vertices[3*triangle_id + ((i + 2)%3)];
i_previous[1] = i_next[0];
}
} /* End of loop on triangle's vertices */
} /* End of loop on triangles on each side of edge */
/* Test quadrilateral's convexity */
if ( _polygon_vertex_is_convex(i_previous[0],
i_0,
i_next[0],
coords) == true
&& _polygon_vertex_is_convex(i_previous[1],
i_1,
i_next[1],
coords) == true)
convex_quad = true;
else
convex_quad = false;
/* Check if edge is locally Delaunay */
edge_locally_delaunay =
_edge_is_locally_delaunay(i_0,
i_1,
vertex_flip[0],
vertex_flip[1],
coords);
/* If edge is not locally Delaunay */
/*---------------------------------*/
if ( edge_locally_delaunay == false
&& convex_quad == true) {
i_min = FVM_MIN(vertex_flip[0], vertex_flip[1]);
i_max = FVM_MAX(vertex_flip[0], vertex_flip[1]);
flip_edge_id = _EDGE_INDEX(i_min, i_max);
for (j = 0; j < 2; j++) { /* Loop on triangles on each side of edge */
triangle_id_0 = edge_neighbors[2*edge_id + j];
triangle_id_1 = edge_neighbors[2*edge_id + ((j + 1)%2)];
/* Redefine adjacent triangles */
for (i = 0; i < 3; i++) {
vertex_id = triangle_vertices[3*triangle_id_0 + i];
if (vertex_id == edge_vertices[2*edge_id + ((j + 1)%2)])
triangle_vertices[3*triangle_id_0 + i] = vertex_flip[(j + 1)%2];
}
/* Redefine triangle so that vertices appear in increasing order */
_triangle_by_sorted_vertices(triangle_vertices + 3*triangle_id_0);
/* Update neighborhood and non locally Delaunay edges */
for (i = 0; i < 3; i++) { /* Loop on triangle's vertices */
i_0 = triangle_vertices[3*triangle_id_0 + i];
i_1 = triangle_vertices[3*triangle_id_0 + ((i + 1)%3)];
i_min = FVM_MIN(i_0, i_1);
i_max = FVM_MAX(i_0, i_1);
current_edge_id = _EDGE_INDEX(i_min, i_max);
if (current_edge_id < edge_id)
restart = true;
/* If current edge is not the new (flip) egde ... */
if (current_edge_id != flip_edge_id) {
if (edge_neighbors[2*current_edge_id] == triangle_id_1)
edge_neighbors[2*current_edge_id] = triangle_id_0;
else if (edge_neighbors[2*current_edge_id + 1] == triangle_id_1)
edge_neighbors[2*current_edge_id + 1] = triangle_id_0;
/* ... and is not either a boundary edge */
if (edge_neighbors[2*current_edge_id + 1] != -1)
edge_is_delaunay[current_edge_id] = false;
}
} /* End of loop on triangle's vertices */
} /* End of loop on triangles on each side of edge */
triangle_id_0 = edge_neighbors[2*edge_id];
triangle_id_1 = edge_neighbors[2*edge_id + 1];
edge_neighbors[2*flip_edge_id] = triangle_id_0;
edge_neighbors[2*flip_edge_id + 1] = triangle_id_1;
edge_neighbors[2*edge_id] = -1;
edge_neighbors[2*edge_id + 1] = -1;
} /* End if edge was not locally Delaunay */
} /* End if edge was initially supposed non locally Delaunay */
if (edge_id == n_edges - 1 && restart == true) {
restart = false;
edge_id = 0;
}
else if (edge_id == n_edges - 1 && restart == false)
face_is_delaunay = true;
else
edge_id++;
} /* End of flip algorithm */
}
/*============================================================================
* Public function definitions
*============================================================================*/
/*----------------------------------------------------------------------------
* Create a structure necessary to the polygon triangulation algorithm.
*
* parameters:
* n_vertices_max <-- maximum expected number of vertices per polygon.
*
* returns:
* pointer to polygon triangulation state structure.
*----------------------------------------------------------------------------*/
fvm_triangulate_state_t *
fvm_triangulate_state_create(const int n_vertices_max)
{
fvm_triangulate_state_t *this_state = NULL;
int n_edges_max = (2*n_vertices_max) - 3;
int n_edges_tot_max = n_edges_max * (n_edges_max - 1) / 2;
BFT_MALLOC(this_state, 1, fvm_triangulate_state_t);
if (n_vertices_max > 3) {
BFT_MALLOC(this_state->triangle_vertices, (n_vertices_max - 2) * 3, int);
BFT_MALLOC(this_state->coords, n_vertices_max*3, fvm_coord_t);
BFT_MALLOC(this_state->list_previous, n_vertices_max, int);
BFT_MALLOC(this_state->list_next, n_vertices_max, int);
BFT_MALLOC(this_state->edge_vertices, n_edges_tot_max*2, int);
BFT_MALLOC(this_state->edge_neighbors, n_edges_tot_max*2, int);
BFT_MALLOC(this_state->edge_is_delaunay, n_edges_tot_max, _Bool);
BFT_MALLOC(this_state->concave, n_vertices_max, _Bool);
}
else {
this_state->triangle_vertices = NULL;
this_state->coords = NULL;
this_state->list_previous = NULL;
this_state->list_next = NULL;
this_state->edge_vertices = NULL;
this_state->edge_neighbors = NULL;
this_state->edge_is_delaunay = NULL;
this_state->concave = NULL;
}
this_state->n_vertices_max = n_vertices_max;
return this_state;
}
/*----------------------------------------------------------------------------
* Destroy a structure necessary to the polygon triangulation algorithm.
*
* parameters:
* this_state <-> pointer to structure that should be destroyed.
*
* returns:
* NULL pointer.
*----------------------------------------------------------------------------*/
fvm_triangulate_state_t *
fvm_triangulate_state_destroy(fvm_triangulate_state_t *this_state)
{
if (this_state != NULL) {
if (this_state->triangle_vertices != NULL) {
BFT_FREE(this_state->triangle_vertices);
BFT_FREE(this_state->coords);
BFT_FREE(this_state->list_previous);
BFT_FREE(this_state->list_next);
BFT_FREE(this_state->edge_vertices);
BFT_FREE(this_state->edge_neighbors);
BFT_FREE(this_state->edge_is_delaunay);
BFT_FREE(this_state->concave);
}
BFT_FREE(this_state);
}
return NULL;
}
/*----------------------------------------------------------------------------
* Triangulate a polygonal face.
*
* For a polygon with n vertices, we should obtain a triangluation with
* (n-2) triangles and (2n-3) edges. If the polygon_vertices argument
* is NULL, 1, 2, ...,n local numbering is implied.
*
* parameters:
* dim <-- spatial dimension (2 or 3).
* n_vertices <-- number of vertices defining the polygon.
* coords <-- coordinates of the triangulation's vertices.
* parent_vertex_num <-- optional indirection to vertex coordinates (1 to n).
* polygon_vertices <-- polygon connectivity; size: n_vertices or empty.
* triangle_vertices --> triangles connectivity;
* size: (n_vertices - 2) * 3.
* state <-> associated triangulation state structure.
*
* returns:
* number of resulting triangles.
*----------------------------------------------------------------------------*/
int
fvm_triangulate_polygon(const int dim,
const int n_vertices,
const fvm_coord_t coords[],
const fvm_lnum_t parent_vertex_num[],
const fvm_lnum_t polygon_vertices[],
fvm_lnum_t triangle_vertices[],
fvm_triangulate_state_t *const state)
{
int i, j;
int n_triangles = 0;
int n_tries = 0;
double epsilon[] = {1.0e-1, 1.0e-2, 0.0, -1.0e-2, -1.0e-1};
int *const list_previous = state->list_previous;
int *const list_next = state->list_next;
_Bool *const concave = state->concave;
/* Initialize state structure */
if (parent_vertex_num != NULL) {
if (polygon_vertices != NULL) {
for (i = 0; i < n_vertices; i++) {
int vertex_id = parent_vertex_num[polygon_vertices[i]-1] - 1;
for (j = 0; j < dim; j++)
state->coords[i*dim + j] = coords[vertex_id*dim + j];
}
}
else {
for (i = 0; i < (n_vertices * dim); i++) {
int vertex_id = parent_vertex_num[i] - 1;
state->coords[i] = coords[vertex_id];
}
}
}
else { /* (if parent_vertex_num == NULL) */
if (polygon_vertices != NULL) {
for (i = 0; i < n_vertices; i++) {
for (j = 0; j < dim; j++)
state->coords[i*dim + j] = coords[(polygon_vertices[i]-1)*dim + j];
}
}
else {
for (i = 0; i < (n_vertices * dim); i++)
state->coords[i] = coords[i];
}
}
/* Determine the work plane (3d coords are overwritten) */
if (dim == 3)
_polygon_plane_3d(n_vertices, state->coords);
/* Initialization */
while ((n_triangles != (n_vertices - 2)) &&
n_tries < 5) {
n_triangles = 0; /* Reset if previous try was a failure */
for (i = 0; i < n_vertices; i++) {
list_previous[i] = i - 1;
list_next[i] = i + 1;
}
list_previous[0] = n_vertices - 1;
list_next[n_vertices - 1] = 0;
for (i = 0; i < n_vertices; i++) {
if (_polygon_vertex_is_convex(list_previous[i],
i,
list_next[i],
state->coords) == true)
concave[i] = false;
else
concave[i] = true;
}
i = 2;
while (i != 0 && i != n_vertices) {
if (_polygon_vertex_is_ear(n_vertices,
list_previous[i],
list_previous,
list_next,
state->concave,
state->coords,
epsilon[n_tries]) == true) {
/* Add a triangle with vertices list_previous[list_previous[i]],
list_previous[i], i */
state->triangle_vertices[n_triangles*3 ] = list_previous
[list_previous[i]];
state->triangle_vertices[n_triangles*3 + 1] = list_previous[i];
state->triangle_vertices[n_triangles*3 + 2] = i;
n_triangles += 1;
/* Cut the ear corresponding to list_previous[i] */
list_previous[i] = list_previous[list_previous[i]];
list_next[list_previous[i]] = i;
if ( (concave[i] == true)
&& (_polygon_vertex_is_convex(list_previous[i],
i,
list_next[i],
state->coords) == true))
concave[i] = false;
if ( (concave[list_previous[i]] == true)
&& (_polygon_vertex_is_convex(list_previous[list_previous[i]],
list_previous[i],
i,
state->coords) == true))
concave[list_previous[i]] = false;
if (list_previous[i] == 0)
i = list_next[i];
}
else /* ! _polygon_vertex_is_ear(...) */
i = list_next[i];
}
n_tries++;
}
/* Now that we have an initial triangulation, apply flip algorithm
to obtain a Delaunay triangulation */
if (n_triangles == n_vertices - 2)
_polygon_delaunay_flip(n_vertices,
state->triangle_vertices,
state->edge_vertices,
state->edge_neighbors,
state->edge_is_delaunay,
state->coords);
/* Update triangle_vertices argument */
if (polygon_vertices != NULL) {
for (i = 0; i < n_triangles * 3; i++)
triangle_vertices[i] = polygon_vertices[state->triangle_vertices[i]];
}
else {
for (i = 0; i < n_triangles * 3; i++)
triangle_vertices[i] = state->triangle_vertices[i] + 1;
}
return n_triangles;
}
/*----------------------------------------------------------------------------*/
#ifdef __cplusplus
}
#endif /* __cplusplus */
syntax highlighted by Code2HTML, v. 0.9.1