Commit 4665e604 authored by Célestin Marot's avatar Célestin Marot
Browse files

update that handles internal triangles better

parent fb1862f2
......@@ -80,9 +80,8 @@ novelties:
can be a totally arbitrary user-defined function. In comparison, Liu et al.
use measure of angles that won't work with anysotropic meshes.
- **The order of candidate nodes** is chosen regarding multiple criteria that
are very different from Liu et al. A candidate node has is chosen in priority
if:
- It is the opposite node in one, two or all adjacent triangles
are very different from Liu et al. A candidate node is chosen in priority if:
- It is the opposite node in 1, 2 or all 3 adjacent triangles
- It is an interior node
- It appears in a lot of triangles
- It results in a good tetrahedron
......
......@@ -252,10 +252,8 @@ static inline int tri_edge_intersection(int t012v0,
static inline int tet_tri_intersection(SPRCavity* SPR,
SPRBbox* tetBbox,
uint8_t tet[4],
unsigned triangleID)
uint8_t* tri)
{
uint8_t* tri = SPR->triangles.array[triangleID].node;
SPRBbox triangleBbox;
SPRBbox_from(&triangleBbox, SPR->points.array[tri[0]].coord);
SPRBbox_add(&triangleBbox, SPR->points.array[tri[1]].coord);
......@@ -463,9 +461,8 @@ static inline int tet_tri_intersection(SPRCavity* SPR,
static SPRNOINLINE int tet_edge_intersection(SPRCavity* SPR,
SPRBbox* tetBbox,
uint8_t tet[4],
unsigned edgeID)
uint8_t* edge)
{
uint8_t* edge = SPR->edges.array[edgeID].node;
double* l0 = SPR->points.array[edge[0]].coord;
double* l1 = SPR->points.array[edge[1]].coord;
if((tetBbox->max[0]<l0[0] && tetBbox->max[0]<l1[0]) ||
......@@ -607,32 +604,6 @@ static SPRNOINLINE double add_quality_map(SPRCavity* SPR,
double* c1 = SPR->points.array[v1].coord;
double* c2 = SPR->points.array[v2].coord;
double* c3 = SPR->points.array[v3].coord;
{
SPRBbox bbox;
SPRBbox_from(&bbox, c0);
SPRBbox_add(&bbox, c1);
SPRBbox_add(&bbox, c2);
SPRBbox_add(&bbox, c3);
uint8_t tet[4] = {v0, v1, v2, v3};
const int n = SPR->points.num;
for (int ptID=0; ptID<n; ptID++) {
if(tet_contains_pt(SPR, &bbox, tet, ptID)) {
SPR->map.qualities[index] = -DBL_MAX;
return -DBL_MAX;
}
}
for (int edgeID=0; edgeID<SPR->edges.num; edgeID++) {
if(tet_edge_intersection(SPR, &bbox, tet, edgeID)) {
SPR->map.qualities[index] = -DBL_MAX;
return -DBL_MAX;
}
}
}
double qual = SPR->quality.function(c0, c1, c2, c3,
SPR->quality.userData);
......@@ -651,22 +622,25 @@ static double get_quality_map(SPRCavity* SPR,
unsigned v0,
unsigned v1,
unsigned v2,
unsigned v3)
unsigned v3,
unsigned* index_ptr)
{
int orientation = get_orient3d(SPR, v0, v1, v2, v3);
if(orientation > 0)
return -1.0;
else if(orientation==0)
return 0.0;
else if(SPR->quality.function==NULL)
return 1.0;
if(orientation >= 0.0)
return -DBL_MAX;
unsigned index = get_compressed_index(v0, v1, v2, v3);
if(index_ptr!=NULL)
*index_ptr = index;
double qual = SPR->map.qualities[index];
if(!isnan(qual)) {
return qual;
}
else if(SPR->quality.function==NULL) {
return DBL_MAX;
}
else {
return add_quality_map(SPR, v0, v1, v2, v3,
index);
......@@ -680,7 +654,7 @@ static double get_quality_map(SPRCavity* SPR,
* that the best tetrahedralization could have */
static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
{
const uint16_t nfaces = SPR->triangles.num;
const uint16_t nfaces = SPR->bndTriangles.num;
const uint8_t npts = SPR->points.num;
int minValid = 255;
......@@ -690,7 +664,7 @@ static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
// new faces have a better chances of being bad, that why we decrement
for (uint16_t i=nfaces; i>0; i--) {
uint16_t triangleID = i-1;
SPRTet tet = SPR->triangles.array[triangleID];
SPRTet tet = SPR->bndTriangles.array[triangleID];
int numValid = 0;
double max = 0.0;
......@@ -707,7 +681,7 @@ static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
tet.node[0],
tet.node[1],
tet.node[2],
tet.node[3]);
tet.node[3], NULL);
if(qual<=SPR->tetrahedra.quality)
continue;
......@@ -757,7 +731,7 @@ static inline void compute_candidates(SPRCavity* SPR,
double qual = get_quality_map(SPR, tet->node[0],
tet->node[1],
tet->node[2],
tet->node[3]);
tet->node[3], NULL);
if(qual<=SPR->tetrahedra.quality)
continue;
......@@ -831,28 +805,28 @@ static void add_face(SPRCavity* SPR,
uint8_t v2)
{
// use compound literal operator (C99 feature)
SPR->triangles.array[SPR->triangles.num].node[0] = v0;
SPR->triangles.array[SPR->triangles.num].node[1] = v1;
SPR->triangles.array[SPR->triangles.num].node[2] = v2;
SPR->triangles.array[SPR->triangles.num].node[3] = v0;
add_face_map(SPR, v0, v1, v2, SPR->triangles.num);
SPR->bndTriangles.array[SPR->bndTriangles.num].node[0] = v0;
SPR->bndTriangles.array[SPR->bndTriangles.num].node[1] = v1;
SPR->bndTriangles.array[SPR->bndTriangles.num].node[2] = v2;
SPR->bndTriangles.array[SPR->bndTriangles.num].node[3] = v0;
add_face_map(SPR, v0, v1, v2, SPR->bndTriangles.num);
increment_3valences(SPR, v0, v1, v2);
SPR->triangles.num++;
SPR->bndTriangles.num++;
}
static void remove_face(SPRCavity* SPR, uint16_t triangleID)
{
SPRTriangle* del = &SPR->triangles.array[triangleID];
SPRTriangle* del = &SPR->bndTriangles.array[triangleID];
remove_face_map(SPR, del->node[0], del->node[1], del->node[2]);
decrement_3valences(SPR, del->node[0], del->node[1], del->node[2]);
if(triangleID!=SPR->triangles.num-1) {
SPRTriangle* last = &SPR->triangles.array[SPR->triangles.num-1];
if(triangleID!=SPR->bndTriangles.num-1) {
SPRTriangle* last = &SPR->bndTriangles.array[SPR->bndTriangles.num-1];
add_face_map(SPR, last->node[0], last->node[1], last->node[2], triangleID);
*del = *last;
}
SPR->triangles.num--;
SPR->bndTriangles.num--;
}
static inline void new_tetrahedralization(SPRCavity* SPR, SPRStep* step)
......@@ -885,12 +859,8 @@ static inline void new_tetrahedralization(SPRCavity* SPR, SPRStep* step)
/***************************
* Initialization
***************************/
static inline void SPR_init(SPRCavity* SPR)
static void SPR_clear_maps(SPRCavity* SPR)
{
assert(SPR->points.num <= SPR_MAX_PTS);
assert(SPR->edges.num <= SPR_MAX_EDGES);
assert(SPR->triangles.num <= SPR_MAX_FACES);
{ // initializing map.face
size_t face_len = sizeof(SPR->map.faces)/sizeof(uint16_t);
......@@ -919,29 +889,44 @@ static inline void SPR_init(SPRCavity* SPR)
}
}
#endif
}
static void SPR_detect_interior_points(SPRCavity* SPR)
{
for (int i=0; i<SPR->points.num; i++) {
if(SPR->points.array[i].valence==0) {
SPR->points.array[i].is_interior = 1;
/* interior points should have a positive valence */
SPR->points.array[i].valence = 1;
}
}
}
static void SPR_init_faceMap_and_valences(SPRCavity* SPR)
{
for (int i=0; i<SPR->points.num; i++) {
SPRPoint* p = &SPR->points.array[i];
p->valence = 0;
p->is_interior = 0;
}
for (uint16_t i=0; i<SPR->triangles.num; i++) {
uint8_t* v = SPR->triangles.array[i].node;
for (uint16_t i=0; i<SPR->bndTriangles.num; i++) {
uint8_t* v = SPR->bndTriangles.array[i].node;
v[3] = v[0];
increment_3valences(SPR, v[0], v[1], v[2]);
add_face_map(SPR, v[0], v[1], v[2], i);
}
}
/* interior points should have a positive valence */
for (int i=0; i<SPR->points.num; i++) {
if(SPR->points.array[i].valence==0) {
SPR->points.array[i].is_interior = 1;
SPR->points.array[i].valence = 1;
}
else {
SPR->points.array[i].is_interior = 0;
}
}
static void SPR_init_step(SPRCavity* SPR)
{
SPR->steps.num = 0;
SPRStep* step = &SPR->steps.array[0];
step->numCandidates = -1; // candidates not computed yet
step->quality = DBL_MAX; // we can hope for the biggest quality at the moment
}
......@@ -965,22 +950,19 @@ static inline void SPR_init(SPRCavity* SPR)
|*| | | \ \ / -_) _/ _ \ ' \| ' \/ -_) _| _| / _ \ ' \ |*|
|*| |_| \_\\___\__\___/_||_|_||_\___\__|\__|_\___/_||_| |*|
\*|=======================================================|*/
int hxtSPR(SPRCavity* SPR)
{
SPR_init(SPR);
SPR->steps.num = 0;
SPRStep* step = &SPR->steps.array[0];
step->numCandidates = -1; // candidates not computed yet
step->quality = DBL_MAX; // we can hope for the biggest quality at the moment
static int hxtSPR_advanced(SPRCavity* SPR) {
// this variable could be avoided, as step is
// SPR->steps.array[SPR->steps.num] all along
SPRStep* step = &SPR->steps.array[SPR->steps.num];
while(1) {
if(SPR->num_search_nodes >= SPR->max_search_nodes){
return 1;
}
int rewind = 0;
if(SPR->triangles.num==0) { // we found a new (better) tetrahedralization
if(SPR->bndTriangles.num==0) { // we found a new (better) tetrahedralization
new_tetrahedralization(SPR, step);
rewind = 1;
}
......@@ -995,7 +977,7 @@ int hxtSPR(SPRCavity* SPR)
rewind = 1;
}
else {
step->tet = SPR->triangles.array[triangleID];
step->tet = SPR->bndTriangles.array[triangleID];
// create (and order) the candidate array
compute_candidates(SPR, step, qualMax);
......@@ -1008,11 +990,12 @@ int hxtSPR(SPRCavity* SPR)
if(step->numCandidates>0) {
step->tet.node[3] = step->candidate[--step->numCandidates];
unsigned index;
double qual = get_quality_map(SPR,
step->tet.node[0],
step->tet.node[1],
step->tet.node[2],
step->tet.node[3]);
step->tet.node[3], &index);
// best quality could have changed if there was a rewind
if(qual <= SPR->tetrahedra.quality)
......@@ -1025,10 +1008,45 @@ int hxtSPR(SPRCavity* SPR)
SPRBbox_add(&tetBbox, SPR->points.array[step->tet.node[2]].coord);
SPRBbox_add(&tetBbox, SPR->points.array[step->tet.node[3]].coord);
for (int ptID=0; ptID<SPR->points.num; ptID++) {
if(tet_contains_pt(SPR, &tetBbox, step->tet.node, ptID)) {
stop = 1;
SPR->map.qualities[index] = -DBL_MAX;
break;
}
}
if(stop)
continue;
for (int edgeID=0; edgeID<SPR->CIEdges.num; edgeID++) {
uint8_t* edge = SPR->CIEdges.array[edgeID].node;
if(tet_edge_intersection(SPR, &tetBbox, step->tet.node, edge)) {
stop = 1;
SPR->map.qualities[index] = -DBL_MAX;
break;
}
}
if(stop)
continue;
for (int triangleID=0; triangleID<SPR->CITriangles.num; triangleID++) {
uint8_t* tri = SPR->CITriangles.array[triangleID].node;
if(tet_tri_intersection(SPR, &tetBbox, step->tet.node, tri)) {
stop = 1;
SPR->map.qualities[index] = -DBL_MAX;
break;
}
}
if(stop)
continue;
// test intersection with boundary triangles
for (int otherFaceID=0; otherFaceID<SPR->triangles.num; otherFaceID++)
{
if(tet_tri_intersection(SPR, &tetBbox, step->tet.node, otherFaceID)){
for (int otherFaceID=0; otherFaceID<SPR->bndTriangles.num; otherFaceID++) {
uint8_t* tri = SPR->bndTriangles.array[otherFaceID].node;
if(tet_tri_intersection(SPR, &tetBbox, step->tet.node, tri)){
stop = 1;
break;
}
......@@ -1097,4 +1115,15 @@ int hxtSPR(SPRCavity* SPR)
}
}
}
}
\ No newline at end of file
}
int hxtSPR(SPRCavity* SPR)
{
SPR_clear_maps(SPR);
SPR_init_faceMap_and_valences(SPR);
SPR_detect_interior_points(SPR);
SPR_init_step(SPR);
return hxtSPR_advanced(SPR);
}
......@@ -17,9 +17,8 @@
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#ifndef _HXTSPR_
#define _HXTSPR_
#ifndef HXTSPR_H
#define HXTSPR_H
#ifdef __cplusplus
extern "C" {
......@@ -39,10 +38,22 @@ extern "C" {
// defines to set/unset the main optimizations
#define SPR_SAVE_ORIENT3D
#define SPR_MAX_PTS 32
#define SPR_MAX_EDGES ((SPR_MAX_PTS*SPR_MAX_PTS - SPR_MAX_PTS)/2)
#define SPR_MAX_TETS ((SPR_MAX_PTS*SPR_MAX_PTS - 3*SPR_MAX_PTS - 2)/2)
#define SPR_MAX_FACES (2*SPR_MAX_TETS+2)
/* From the book of Cheng et al. 'Delaunay Mesh Generation'
* A triangulation of n vertices in R^3 has at most `ntet=(n^2 − 3n − 2)/2` tetrahedra,
* at most nedg=(n^2 − n)/2` edges (this is simply the 2-combination of n)
* and at most `ntri=n^2 − 3n = 2*ntet+2` triangles */
#define SPR_MAX_TETS ((SPR_MAX_PTS*SPR_MAX_PTS - 3*SPR_MAX_PTS - 2)/2) //463
/* In pratices however, the number of triangles will certainly never be greater than 500
* because the different constraints of the SPR make the upper bound impossible to be reached
* Same argument for constrained interior triangles and edges */
#define SPR_MAX_BNDTRIANGLES 500
#define SPR_MAX_CITRIANGLES 100
#define SPR_MAX_CIEDGES 100
/* structure for the points (32 bytes)*/
......@@ -78,8 +89,6 @@ typedef struct {
double quality; // quality before inserting this tet
} SPRStep;
typedef struct SPRExact_s SPRExact; // an opaque structure...
/* This is the main structure defining your cavity.
* You might want to allocate this heavy structure on the heap ! */
......@@ -98,14 +107,21 @@ typedef struct {
int num;
} points;
/* Constrained edges
/* Constrained Interior edges
* hxtSPR() might shuffle what you put in here */
struct {
SPREdge array[SPR_MAX_EDGES];
SPREdge array[SPR_MAX_CIEDGES];
int num;
} edges;
} CIEdges;
/* Boundary triangles + constrained triangles (appearing twice)
/* Constrained Interior triangles
* hxtSPR() might shuffle what you put in here */
struct {
SPRTriangle array[SPR_MAX_CITRIANGLES];
int num;
} CITriangles;
/* Boundary triangles
*
* Nodes of boundary triangles must appear in a counter-clockwise order
* when seen from the interior of the cavity.
......@@ -113,14 +129,11 @@ typedef struct {
* Actually, there might be multiple closed surfaces but it is better to
* tetrahedralize them separately...
*
* Interior, constrained triangles must appear in both counter-clockwise
* and clockwise order (when seen from anywhere obviously...)
*
* hxtSPR() might shuffle what you put in here */
* hxtSPR() may shuffle what you put in here */
struct {
SPRTriangle array[SPR_MAX_FACES];
SPRTriangle array[SPR_MAX_BNDTRIANGLES];
int num;
} triangles;
} bndTriangles;
/* After the SPR, `tetrahedra.array` will contain the `tetrahedra.num`
* tetrahedra of the best tetrahedralization found, or be left unchanged
......@@ -137,7 +150,7 @@ typedef struct {
double quality; // known lower bound on quality or 0
} tetrahedra;
// quality function and lower bound on quality
// quality function
struct {
double (*function)(double* p0, double* p1, double* p2, double* p3,
void* userData);
......@@ -174,7 +187,7 @@ typedef struct {
* Strategy for 3-D Mesh Generation and Mesh Improvement"
* by Liu Jianfei, Sun Shuli and Wang Dachuan
*
* It was however heavily optimized, as explained in the paper for the
* It was however slightly optimized, as explained in the paper for the
* 28th IMR called "Reviving the Search for Optimal Tetrahedralizations"
* by Marot Célestin, Verhetsel Kilian and Jean-François Remacle.
*
......@@ -189,7 +202,7 @@ typedef struct {
*
* If you only want to find high-quality tetrahedralization, or if you don't
* want to waste time for poor quality result anyway, limit the number of nodes
* in the search tree by setting max_search_nodes somewhere around SPR_MAX_FACES
* in the search tree by setting max_search_nodes somewhere around 500
*
* To call hxtSPR, fill the cavity structure as described in the structure
* declaration, and, if there is no error (return value is 0),
......
......@@ -21,7 +21,7 @@ Author: Célestin Marot (celestin.marot@uclouvain.be) */
* Same benchmark as in the paper *
* "REVIVING THE SEARCH FOR OPTIMAL TETRAHEDRALIZATIONS", *
* section 2.5 *
* + benchmark in Liu's paper: *
* + ~ benchmark in Liu's paper: *
**********************************************************/
#include "HXTSPR.h"
......@@ -95,15 +95,17 @@ double aspectRatio(double a[3], double b[3], double c[3], double d[3],
}
static int benchmark(double coord[][3], int numCoord, uint8_t face[][3], int numFaces, double original_best) {
SPRCavity SPR = {0}; // fully on the stack ^^
static int benchmark(double coord[][3], int numPoints, uint8_t face[][3], int numFaces, double original_best) {
SPRCavity SPR = {0}; // fully on the stack
SPR.quality.function = aspectRatio;
SPR.max_search_nodes = UINT64_MAX; // no limit :p
SPR.CIEdges.num = 0;
SPR.CITriangles.num = 0;
double total_time = 0.0;
uint64_t total_nodes = 0;
for (int i=0; i<numCoord; i++)
for (int i=0; i<numPoints; i++)
{
for (int j=0; j<3; j++)
{
......@@ -113,14 +115,14 @@ static int benchmark(double coord[][3], int numCoord, uint8_t face[][3], int num
for (int i=0; i<numFaces; i++)
{
SPR.triangles.array[i] = (SPRTriangle) {.node = {face[i][0], face[i][1], face[i][2], face[i][0]}};
SPR.bndTriangles.array[i] = (SPRTriangle) {.node = {face[i][0], face[i][1], face[i][2], face[i][0]}};
}
int iter = 0;
for (; iter<100 && total_time < 60.0; iter++) {
SPR.tetrahedra.quality = original_best;
SPR.points.num = numCoord;
SPR.triangles.num = numFaces;
SPR.points.num = numPoints;
SPR.bndTriangles.num = numFaces;
SPR.num_search_nodes = 0;
struct timespec time0, time1;
......@@ -131,9 +133,9 @@ static int benchmark(double coord[][3], int numCoord, uint8_t face[][3], int num
total_nodes += SPR.num_search_nodes;
#ifndef NDEBUG
if(SPR.triangles.num!=numFaces) {
if(SPR.bndTriangles.num!=numFaces) {
fprintf(stderr, "%d triangles before, %d after\n",
numFaces, SPR.triangles.num);
numFaces, SPR.bndTriangles.num);
return 4;
}
......@@ -141,18 +143,18 @@ static int benchmark(double coord[][3], int numCoord, uint8_t face[][3], int num
for (int j=0; j<numFaces; j++) {
int found = 0;
for (int i=0; i<numFaces; i++) {
if((SPR.triangles.array[i].node[0]==face[j][0] &&
SPR.triangles.array[i].node[1]==face[j][1] &&
SPR.triangles.array[i].node[2]==face[j][2] &&
SPR.triangles.array[i].node[3]==SPR.triangles.array[i].node[0]) ||
(SPR.triangles.array[i].node[0]==face[j][1] &&
SPR.triangles.array[i].node[1]==face[j][2] &&
SPR.triangles.array[i].node[2]==face[j][0] &&
SPR.triangles.array[i].node[3]==SPR.triangles.array[i].node[0]) ||
(SPR.triangles.array[i].node[0]==face[j][2] &&
SPR.triangles.array[i].node[1]==face[j][0] &&
SPR.triangles.array[i].node[2]==face[j][1] &&
SPR.triangles.array[i].node[3]==SPR.triangles.array[i].node[0]) )
if((SPR.bndTriangles.array[i].node[0]==face[j][0] &&
SPR.bndTriangles.array[i].node[1]==face[j][1] &&
SPR.bndTriangles.array[i].node[2]==face[j][2] &&
SPR.bndTriangles.array[i].node[3]==SPR.bndTriangles.array[i].node[0]) ||
(SPR.bndTriangles.array[i].node[0]==face[j][1] &&
SPR.bndTriangles.array[i].node[1]==face[j][2] &&
SPR.bndTriangles.array[i].node[2]==face[j][0] &&
SPR.bndTriangles.array[i].node[3]==SPR.bndTriangles.array[i].node[0]) ||
(SPR.bndTriangles.array[i].node[0]==face[j][2] &&
SPR.bndTriangles.array[i].node[1]==face[j][0] &&
SPR.bndTriangles.array[i].node[2]==face[j][1] &&
SPR.bndTriangles.array[i].node[3]==SPR.bndTriangles.array[i].node[0]) )
{
found = 1;
break;
......@@ -697,7 +699,7 @@ int main()
uint8_t (* face)[3] = malloc(3*sizeof(uint8_t) * (12+(32-8)*2));
// we begin with a cube
int numCoord = 8;
int numPoints = 8;
double c = 1/sqrt(3);
memcpy(coord[0], (double[3]){-c,-c,-c}, 3*sizeof(double) );
memcpy(coord[1], (double[3]){ c,-c,-c}, 3*sizeof(double) );
......@@ -724,8 +726,8 @@ int main()
printf("\nusing growing number of points\n");
for (int i=8; i<=32; i++) {
printf("with %d points, %d faces\n", numCoord, numFaces);
HXT_CHECK( benchmark(coord, numCoord, face, numFaces, 0.0) );
printf("with %d points, %d faces\n", numPoints, numFaces);
HXT_CHECK( benchmark(coord, numPoints, face, numFaces, 0.0) );
if(i==32)
continue;
......@@ -759,7 +761,7 @@ int main()
uint8_t v0 = face[biggestTri][0];
uint8_t v1 = face[biggestTri][1];
uint8_t v2 = face[biggestTri][2];
uint8_t v3 = numCoord++;
uint8_t v3 = numPoints++;
// we add a point on the sphere in the middle of the triangle
double *v0Coord = coord[v0];
......
......@@ -210,12 +210,13 @@ int main()
SPR->points.array[i].coord[j] = coord[i][j];
}
SPR->triangles.num = sizeof(face)/sizeof(face[0]);
for (int i=0; i<SPR->triangles.num; i++) {
SPR->triangles.array[i] = (SPRTriangle) {.node = {face[i][0], face[i][1], face[i][2]}};
SPR->bndTriangles.num = sizeof(face)/sizeof(face[0]);
for (int i=0; i<SPR->bndTriangles.num; i++) {
SPR->bndTriangles.array[i] = (SPRTriangle) {.node = {face[i][0], face[i][1], face[i][2]}};
}