Commit 6e889870 authored by Célestin Marot's avatar Célestin Marot
Browse files

now the benchmark include a growing polyhedron

parent 79ad3557
......@@ -26,6 +26,7 @@ Author: Célestin Marot (celestin.marot@uclouvain.be) */
#include "HXTSPR.h"
#include <math.h>
#include <time.h>
#include <string.h>
#define HXT_CHECK(x) do{ int _err=x; if(_err) { fprintf(stderr, "error %d\n",_err); return _err;}}while(0)
......@@ -92,12 +93,13 @@ double aspectRatio(double a[3], double b[3], double c[3], double d[3],
}
static int benchmark(double coord[][3], int numCoord, unsigned char face[][3], int numFaces, double original_best) {
static int benchmark(double coord[][3], int numCoord, 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
double total_time = 0.0;
uint64_t total_nodes = 0;
for (int i=0; i<numCoord; i++)
{
......@@ -124,6 +126,7 @@ static int benchmark(double coord[][3], int numCoord, unsigned char face[][3], i
HXT_CHECK( hxtSPR(&SPR) );
timespec_get(&time1, TIME_UTC);
total_time += (time1.tv_sec - time0.tv_sec) + 0.000000001*(time1.tv_nsec - time0.tv_nsec);
total_nodes += SPR.num_search_nodes;
#ifndef NDEBUG
if(SPR.triangles.num!=numFaces) {
......@@ -166,8 +169,8 @@ static int benchmark(double coord[][3], int numCoord, unsigned char face[][3], i
fflush(stdout);
}
printf("%.9f sec in average, from quality %f to %f\n",
total_time/iter, original_best, SPR.tetrahedra.quality);
printf("%.9f sec, %.2f search nodes, quality %f -> %f\n",
total_time/iter, (double) total_nodes/iter, original_best, SPR.tetrahedra.quality);
return 0;
}
......@@ -204,7 +207,7 @@ int main()
{0.59999999999999975575093458246556110680103302001953,-0.00000000000000014695761589768230597476116143382705,0.17320508075688789695334435236873105168342590332031}
};
unsigned char face[44][3] = {
uint8_t face[44][3] = {
{2, 1, 4},
{4, 1, 6},
{3, 4, 6},
......@@ -294,7 +297,7 @@ int main()
};
unsigned char face[52][3] = {
uint8_t face[52][3] = {
{4, 6, 3},
{4, 3, 9},
{5, 10, 3},
......@@ -390,7 +393,7 @@ int main()
{-0.22660211555983006403991453225899022072553634643555,0.15670939947857920326867997573572210967540740966797,0.00244676225018409555789222764587975689209997653961}
};
unsigned char face[52][3] = {
uint8_t face[52][3] = {
{3, 5, 4},
{9, 3, 4},
{5, 3, 11},
......@@ -479,7 +482,7 @@ int main()
{-0.15034402621635831187596465952083235606551170349121,-0.45463342746132362481858990577165968716144561767578,0.14389273927417040033205353211087640374898910522461},
};
unsigned char face[38][3] = {
uint8_t face[38][3] = {
{1, 4, 2},
{1, 2, 7},
{2, 4, 8},
......@@ -560,7 +563,7 @@ int main()
{0.34251857348020220417339487539720721542835235595703,0.80587333785436421340619972397689707577228546142578,0.96361204438398462279735667834756895899772644042969}
};
unsigned char face[48][3] = {
uint8_t face[48][3] = {
{8, 3, 4},
{9, 4, 3},
{3, 5, 11},
......@@ -627,7 +630,7 @@ int main()
{0.8,0.8,0.8},
};
unsigned char face[6][3] = {
uint8_t face[6][3] = {
{0, 1, 2},
{0, 2, 3},
{0, 3, 1},
......@@ -638,7 +641,7 @@ int main()
double best = aspectRatio(coord[1], coord[2], coord[3], coord[4], NULL);
printf("\nused for 2-3 flip\n");
printf("\nused for 2-3 or 3-2 flip\n");
HXT_CHECK( benchmark(coord, sizeof(coord)/sizeof(coord[0]), face, sizeof(face)/sizeof(face[0]), best) );
HXT_CHECK( benchmark(coord, sizeof(coord)/sizeof(coord[0]), face, sizeof(face)/sizeof(face[0]), 0.0) );
}
......@@ -656,7 +659,7 @@ int main()
{0, -0.1, 1.0}
};
unsigned char face[12][3] = {
uint8_t face[12][3] = {
{0, 1, 6},
{1, 2, 6},
{2, 3, 6},
......@@ -671,7 +674,7 @@ int main()
{5, 7, 0},
};
// I don't remember the worst one
// I don't remember the worst one...
double best = 1000.0;
for(int i=0; i<6; i++) {
double qual = aspectRatio(coord[i], coord[(i+1)%6], coord[6], coord[7], NULL);
......@@ -684,5 +687,100 @@ int main()
HXT_CHECK( benchmark(coord, sizeof(coord)/sizeof(coord[0]), face, sizeof(face)/sizeof(face[0]), 0.0) );
}
{
// for different size of polyhedron
double (* coord)[3] = malloc(3*sizeof(double) * 40); // we will go up to 40 points
uint8_t (* face)[3] = malloc(3*sizeof(uint8_t) * (12+(40-8)*2));
// we begin with a cube
int numCoord = 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) );
memcpy(coord[2], (double[3]){ c, c,-c}, 3*sizeof(double) );
memcpy(coord[3], (double[3]){-c, c,-c}, 3*sizeof(double) );
memcpy(coord[4], (double[3]){-c,-c, c}, 3*sizeof(double) );
memcpy(coord[5], (double[3]){ c,-c, c}, 3*sizeof(double) );
memcpy(coord[6], (double[3]){ c, c, c}, 3*sizeof(double) );
memcpy(coord[7], (double[3]){-c, c, c}, 3*sizeof(double) );
int numFaces = 12;
memcpy(face[0], (uint8_t[3]){0,1,2}, 3*sizeof(uint8_t));
memcpy(face[1], (uint8_t[3]){0,2,3}, 3*sizeof(uint8_t));
memcpy(face[2], (uint8_t[3]){7,5,4}, 3*sizeof(uint8_t));
memcpy(face[3], (uint8_t[3]){7,6,5}, 3*sizeof(uint8_t));
memcpy(face[4], (uint8_t[3]){0,4,5}, 3*sizeof(uint8_t));
memcpy(face[5], (uint8_t[3]){0,5,1}, 3*sizeof(uint8_t));
memcpy(face[6], (uint8_t[3]){2,6,7}, 3*sizeof(uint8_t));
memcpy(face[7], (uint8_t[3]){2,7,3}, 3*sizeof(uint8_t));
memcpy(face[8], (uint8_t[3]){7,0,3}, 3*sizeof(uint8_t));
memcpy(face[9], (uint8_t[3]){7,4,0}, 3*sizeof(uint8_t));
memcpy(face[10], (uint8_t[3]){2,5,6}, 3*sizeof(uint8_t));
memcpy(face[11], (uint8_t[3]){2,1,5}, 3*sizeof(uint8_t));
printf("\nusing growing number of points\n");
for (int i=8; i<=40; i++) {
printf("with %d points, %d faces\n", numCoord, numFaces);
HXT_CHECK( benchmark(coord, numCoord, face, numFaces, 0.0) );
if(i==40)
continue;
// we select the biggest triangle and split it by inserting a point
// in its barycenter
int biggestTri = 0;
double maxArea2 = 0.0;
for (int tri=0; tri<numFaces; tri++) {
double *v0Coord = coord[face[tri][0]];
double *v1Coord = coord[face[tri][1]];
double *v2Coord = coord[face[tri][2]];
double ab[3], ac[3], area2=0.0;
for (int j=0; j<3; j++) {
ab[j] = v1Coord[j] - v0Coord[j];
ac[j] = v2Coord[j] - v0Coord[j];
}
for (int j=0; j<3; j++) {
double nj = ab[(j+1)%3]*ac[(j+2)%3] - ab[(j+2)%3]*ac[(j+1)%3];
area2 += nj*nj;
}
if(area2 > maxArea2) {
maxArea2 = area2;
biggestTri = tri;
}
}
uint8_t v0 = face[biggestTri][0];
uint8_t v1 = face[biggestTri][1];
uint8_t v2 = face[biggestTri][2];
uint8_t v3 = numCoord++;
// we add a point on the sphere in the middle of the triangle
double *v0Coord = coord[v0];
double *v1Coord = coord[v1];
double *v2Coord = coord[v2];
double v3Coord[3];
double normalize = 0.0;
for (int j=0; j<3; j++) {
v3Coord[j] = (v0Coord[j] + v1Coord[j] + v2Coord[j])/3.0;
normalize += v3Coord[j]*v3Coord[j];
}
normalize = sqrt(normalize);
for (int j=0; j<3; j++) {
v3Coord[j]/=normalize;
}
memcpy(coord[v3], v3Coord, 3*sizeof(double) );
memcpy(face[biggestTri], (uint8_t[3]){v0,v1,v3}, 3*sizeof(uint8_t));
memcpy(face[numFaces++], (uint8_t[3]){v1,v2,v3}, 3*sizeof(uint8_t));
memcpy(face[numFaces++], (uint8_t[3]){v2,v0,v3}, 3*sizeof(uint8_t));
}
}
return 0;
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment