Commit 82bb3a0b authored by Célestin Marot's avatar Célestin Marot
Browse files

recursive version works, iterative has a bug

parent 62b5fa26
CFLAGS=-Wall -O3 -g -DDEBUG
LDFLAGS=-lm
CFLAGS=-Wall -O3 -g -DDEBUG -Wextra -fno-omit-frame-pointer -Wno-missing-field-initializers -Wpointer-arith -Wcast-align -Wcast-qual -Wstrict-aliasing -Wpointer-arith -Winit-self -Wredundant-decls -fsanitize=address -fsanitize=leak -fsanitize=undefined
LDFLAGS=-lm -fsanitize=address -fsanitize=leak -fsanitize=undefined
EXEC=./exec
SRC= $(wildcard *.c)
......
......@@ -31,7 +31,7 @@ Author: Célestin Marot (celestin.marot@uclouvain.be) */
#endif
// get the index in a tet. of node i from facet j by doing nodeFromFacet[i][j];
const static int nodeFromFacet[3][4] = {{1, 2, 3, 0 },
static const int nodeFromFacet[3][4] = {{1, 2, 3, 0 },
{3, 3, 1, 1 },
{2, 0, 0, 2 }};
......@@ -176,14 +176,14 @@ static SPRNOINLINE uint32_t add_chirotope(SPRCavity* SPR,
ID = v0*n3+v1*n2+v2*n1+v3;
index = ID/16;
bit = ID%16 * 2;
SPR->map.chirotope[index] &= ~(3U<<bit); // clear bits
SPR->map.chirotope[index] |= val<<bit; // set bits
SPR->map.chirotopes[index] &= ~(3U<<bit); // clear bits
SPR->map.chirotopes[index] |= val<<bit; // set bits
ID = v0*n3+v1*n2+v3*n1+v2;
index = ID/16;
bit = ID%16 * 2;
SPR->map.chirotope[index] &= ~(3U<<bit); // clear bits
SPR->map.chirotope[index] |= oppval<<bit; // set bits
SPR->map.chirotopes[index] &= ~(3U<<bit); // clear bits
SPR->map.chirotopes[index] |= oppval<<bit; // set bits
}
return val;
......@@ -201,7 +201,7 @@ static int get_chirotope(SPRCavity* SPR,
unsigned ID = v0*(n*n*n)+v1*(n*n)+v2*n+v3;
unsigned index = ID/16;
unsigned bit = ID%16 * 2;
uint32_t val = SPR->map.chirotope[index]>>bit & 3;
uint32_t val = SPR->map.chirotopes[index]>>bit & 3;
if(val==0) {
return (int)add_chirotope(SPR, v0, v1, v2 ,v3) - 2;
}
......@@ -625,16 +625,16 @@ static SPRNOINLINE double add_qualityMap(SPRCavity* SPR,
SPRBboxAddOne(&bbox, c3);
uint8_t tet[4] = {v0, v1, v2, v3};
const unsigned n = SPR->points.num;
const int n = SPR->points.num;
for (unsigned ptID=0; ptID<n; ptID++) {
for (int ptID=0; ptID<n; ptID++) {
if(tet_contains_pt(SPR, &bbox, tet, ptID)) {
SPR->map.qualities[index] = 0.0;
return 0.0;
}
}
for (unsigned edgeID=0; edgeID<SPR->edges.num; edgeID++) {
for (int edgeID=0; edgeID<SPR->edges.num; edgeID++) {
if(tet_edge_intersection(SPR, &bbox, tet, edgeID)) {
SPR->map.qualities[index] = 0.0;
return 0.0;
......@@ -743,7 +743,7 @@ static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
#ifdef SPR_USE_FACE_HEURISTIC
const uint8_t npts = SPR->points.num;
int worstNumValid = SPR_MAX_PTS;
int worstNumValid = 255;
double worstQualSum = DBL_MAX;
uint16_t worstFaceID = UINT16_MAX;
double minimax = DBL_MAX;
......@@ -757,7 +757,7 @@ static inline uint16_t best_face_heuristic(SPRCavity* SPR, double* qualMax)
double qualSum = 0.0;
double max = 0.0;
for (tet.node[3]=0; tet.node[3] < npts && numValid < worstNumValid; tet.node[3]++)
for (tet.node[3]=0; tet.node[3] < npts && numValid <= worstNumValid; tet.node[3]++)
{
if(tet.node[3] == tet.node[0] ||
tet.node[3] == tet.node[1] ||
......@@ -929,6 +929,8 @@ static inline void compute_candidates(SPRCavity* SPR,
sum = tsum;
}
assert(sum==num);
for (uint8_t i=0; i<num; i++) {
step->candidate[count[score[i]]++] = candidate[i];
}
......@@ -1016,11 +1018,12 @@ static inline int SPR_init(SPRCavity* SPR)
}
{ // allocating & initializing map.face
SPR->map.faces = malloc(sizeof(uint16_t) * n*n*n);
size_t face_len = SPR_MAX_PTS*SPR_MAX_PTS*SPR_MAX_PTS;
SPR->map.faces = malloc(sizeof(uint16_t) * face_len);
if(SPR->map.faces==NULL)
goto allocation_error;
for (size_t i=0; i<n*n*n; i++) {
for (size_t i=0; i<face_len; i++) {
SPR->map.faces[i] = UINT16_MAX;
}
}
......@@ -1029,6 +1032,8 @@ static inline int SPR_init(SPRCavity* SPR)
{ // allocating & initializing map.qualities
size_t qual_len = n*(n-1)*(n-2)*(n-3)/24;
SPR->map.qualities = malloc(sizeof(double)*qual_len);
if(SPR->map.qualities==NULL)
goto allocation_error;
for (size_t i=0; i<qual_len; i++) {
SPR->map.qualities[i] = NAN;
}
......@@ -1038,7 +1043,9 @@ static inline int SPR_init(SPRCavity* SPR)
#ifdef SPR_SAVE_CHIROTOPES
{ // allocating & setting to zero map.chirotopes
size_t chir_len = n*n*n*n/16+1;
SPR->map.chirotope = calloc(chir_len, sizeof(uint32_t));
SPR->map.chirotopes = calloc(chir_len, sizeof(uint32_t));
if(SPR->map.chirotopes==NULL)
goto allocation_error;
}
#endif
......@@ -1091,7 +1098,7 @@ static inline void SPR_terminate(SPRCavity* SPR)
free(SPR->map.qualities);
#endif
#ifdef SPR_SAVE_CHIROTOPES
free(SPR->map.chirotopes);
free(SPR->map.chirotopess);
#endif
}
......@@ -1116,89 +1123,168 @@ static inline void SPR_terminate(SPRCavity* SPR)
|*| | | \ \ / -_) _/ _ \ ' \| ' \/ -_) _| _| / _ \ ' \ |*|
|*| |_| \_\\___\__\___/_||_|_||_\___\__|\__|_\___/_||_| |*|
\*|=======================================================|*/
int hxtSPR(SPRCavity* SPR) {
if(SPR_init(SPR))
return 2;
// int hxtSPR(SPRCavity* SPR) {
// if(SPR_init(SPR))
// return 2;
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
while(1){
if(SPR->num_search_nodes >= SPR->max_search_nodes){
SPR_terminate(SPR);
return 1;
}
SPR->num_search_nodes++;
int rewind = 0;
if(SPR->triangles.num==0) { // we found a new (better) tetrahedralization
new_tetrahedralization(SPR, step);
rewind = 1;
}
else if(step->numCandidates<0) {
assert(step->quality > SPR->tetrahedra.quality);
// 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
// we still have to compute the face and candidates
double qualMax;
const uint16_t triangleID = best_face_heuristic(SPR, &qualMax);
// while(1){
// if(SPR->num_search_nodes >= SPR->max_search_nodes){
// SPR_terminate(SPR);
// return 1;
// }
// SPR->num_search_nodes++;
if(triangleID != UINT16_MAX) {
step->tet = SPR->triangles.array[triangleID];
// create the candidate array
compute_candidates(SPR, step, qualMax, triangleID);
// assert(step==&SPR->steps.array[SPR->steps.num]);
// remove the face from the map
remove_face(SPR, triangleID);
}
else {
rewind = 1;
}
}
else if(step->quality <= SPR->tetrahedra.quality) {
rewind = 1;
}
else if(step->numCandidates==0) {
add_face(SPR,
step->tet.node[0],
step->tet.node[1],
step->tet.node[2]);
rewind = 1;
}
// printf("step %d, %d candidates, quality %f\n", SPR->steps.num, step->numCandidates, step->quality);
// int rewind = 0;
// if(SPR->triangles.num==0) { // we found a new (better) tetrahedralization
// printf("tetrahedralization found ! quality %f\n", step->quality);
// new_tetrahedralization(SPR, step);
// rewind = 1;
// }
// else if(step->quality <= SPR->tetrahedra.quality) {
// rewind = 1;
// }
// else if(step->numCandidates<0) {
// assert(step->quality > SPR->tetrahedra.quality);
// // we still have to compute the face and candidates
// double qualMax;
// const uint16_t triangleID = best_face_heuristic(SPR, &qualMax);
// if(triangleID != UINT16_MAX) {
// step->tet = SPR->triangles.array[triangleID];
// // create the candidate array
// compute_candidates(SPR, step, qualMax, triangleID);
// // there still can be 0 candidates because intersection are only computed here above
// }
// if(step->numCandidates>0) {
// // remove the face from the map
// remove_face(SPR, triangleID);
// }
// else {
// rewind = 1;
// }
// printf("found %d candidates\n", step->numCandidates);
// }
// else if(step->numCandidates==0) {
// add_face(SPR,
// step->tet.node[0],
// step->tet.node[1],
// step->tet.node[2]);
// rewind = 1;
// }
// if(rewind) {
// printf("rewinding\n");
// if(SPR->steps.num==0) {
// SPR_terminate(SPR);
// return 0;
// }
// step = &SPR->steps.array[--SPR->steps.num];
// for (int f=0; f<3; f++) {
// uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
// uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
// uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
// uint16_t index = SPR_get_faceMap(SPR, p0, p2, p1);
// if(index==UINT16_MAX) {
// // add the face to the map
// add_face(SPR, p0, p1, p2);
// }
// else {
// // remove the face from the map
// remove_face(SPR, index);
// }
// }
// continue;
// }
// assert(step->numCandidates>0);
// step->tet.node[3] = step->candidate[--step->numCandidates];
// double qual = get_quality_map(SPR,
// step->tet.node[0],
// step->tet.node[1],
// step->tet.node[2],
// step->tet.node[3]);
// // best.quality could have changed...
// if(qual>SPR->tetrahedra.quality) {
// for (int f=0; f<3; f++) {
// uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
// uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
// uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
// uint16_t index = SPR_get_faceMap(SPR, p0, p1, p2);
// if(index==UINT16_MAX) {
// // add the face to the map
// add_face(SPR, p0, p2, p1);
// }
// else {
// // remove the face from the map
// remove_face(SPR, index);
// }
// }
// SPRStep* nextStep = &SPR->steps.array[++SPR->steps.num];
// nextStep->quality = qual < step->quality ? qual : step->quality;
// nextStep->numCandidates = -1;
// step = nextStep;
// }
// }
// }
void hxtSPR_recur(SPRCavity* SPR) {
if(SPR->num_search_nodes >= SPR->max_search_nodes){
return;
}
SPR->num_search_nodes++;
SPRStep* step = &SPR->steps.array[SPR->steps.num];
if(SPR->triangles.num==0) { // we found a new (better) tetrahedralization
new_tetrahedralization(SPR, step);
return;
}
if(rewind) {
{
double qualMax;
const uint16_t triangleID = best_face_heuristic(SPR, &qualMax);
if(SPR->steps.num==0) {
SPR_terminate(SPR);
return 0;
}
step = &SPR->steps.array[--SPR->steps.num];
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
uint16_t index = SPR_get_faceMap(SPR, p0, p2, p1);
if(index==UINT16_MAX) {
// add the face to the map
add_face(SPR, p0, p1, p2);
}
else {
// remove the face from the map
remove_face(SPR, index);
}
}
continue;
}
if(triangleID == UINT16_MAX)
return;
step->tet = SPR->triangles.array[triangleID];
assert(step->numCandidates>0);
step->tet.node[3] = step->candidate[--step->numCandidates];
// create the candidate array
compute_candidates(SPR, step, qualMax, triangleID);
// remove the face from the map
remove_face(SPR, triangleID);
}
for (; step->numCandidates; step->numCandidates--)
{
step->tet.node[3] = step->candidate[step->numCandidates-1];
double qual = get_quality_map(SPR,
step->tet.node[0],
......@@ -1206,31 +1292,78 @@ int hxtSPR(SPRCavity* SPR) {
step->tet.node[2],
step->tet.node[3]);
// best.quality could have changed...
if(qual>SPR->tetrahedra.quality) {
if(qual <= SPR->tetrahedra.quality)
continue;
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
uint16_t index = SPR_get_faceMap(SPR, p0, p1, p2);
uint16_t index = SPR_get_faceMap(SPR, p0, p1, p2);
if(index==UINT16_MAX) {
// add the face to the map
add_face(SPR, p0, p2, p1);
}
else {
// remove the face from the map
remove_face(SPR, index);
}
if(index==UINT16_MAX) {
// add the face to the map
add_face(SPR, p0, p2, p1);
}
else {
// remove the face from the map
remove_face(SPR, index);
}
}
SPRStep* nextStep = &SPR->steps.array[SPR->steps.num+1];
nextStep->quality = qual < step->quality ? qual : step->quality;
nextStep->numCandidates = -1;
SPR->steps.num++;
hxtSPR_recur(SPR);
SPRStep* nextStep = &SPR->steps.array[++SPR->steps.num];
nextStep->quality = qual < step->quality ? qual : step->quality;
nextStep->numCandidates = -1;
SPR->steps.num--;
step = nextStep;
for (int f=0; f<3; f++) {
uint8_t p0 = step->tet.node[nodeFromFacet[0][f]];
uint8_t p1 = step->tet.node[nodeFromFacet[1][f]];
uint8_t p2 = step->tet.node[nodeFromFacet[2][f]];
uint16_t index = SPR_get_faceMap(SPR, p0, p2, p1);
if(index==UINT16_MAX) {
// add the face to the map
add_face(SPR, p0, p1, p2);
}
else {
// remove the face from the map
remove_face(SPR, index);
}
}
if(step->quality <= SPR->tetrahedra.quality ||
SPR->num_search_nodes >= SPR->max_search_nodes)
break;
}
add_face(SPR,
step->tet.node[0],
step->tet.node[1],
step->tet.node[2]);
}
int hxtSPR(SPRCavity* SPR) {
if(SPR_init(SPR))
return 2;
SPR->steps.num = 0;
SPR->steps.array[0].quality = DBL_MAX; // we can hope for the biggest quality at the moment
hxtSPR_recur(SPR);
SPR_terminate(SPR);
if(SPR->num_search_nodes>=SPR->max_search_nodes)
return 1;
return 0;
}
\ No newline at end of file
......@@ -35,7 +35,7 @@ extern "C" {
#define SPR_SAVE_QUALITIES
#define SPR_USE_FACE_HEURISTIC
// #define SPR_SAVE_CHIROTOPES
// #define SPR_INPUT_VERIFICATION
// #define SPR_INPUT_VERIFICATION // not implemented yet
#define SPR_MAX_PTS 48 // you can change only this value
#define SPR_MAX_EDGES ((SPR_MAX_PTS*SPR_MAX_PTS - SPR_MAX_PTS)/2)
......@@ -170,7 +170,7 @@ typedef struct {
double* qualities;
#endif
#ifdef SPR_SAVE_CHIROTOPES
uint32_t* chirotope;
uint32_t* chirotopes;
#endif
uint16_t* faces;
} map;
......
......@@ -120,7 +120,7 @@ int gmshTetDraw(SPRCavity* SPR, const char* filename)
int main(int argc, char const *argv[])
int main()
{
double coord[25][3] = {
{0.58671240234164312443709832223248668015003204345703,0.03516283877473333485674800158449215814471244812012,0.04997715958560345478334241420270700473338365554810},
......@@ -198,6 +198,8 @@ int main(int argc, char const *argv[])
};
SPRCavity* SPR = malloc(sizeof(SPRCavity));
if(SPR==NULL)
return 1;
SPR->max_search_nodes = UINT64_MAX;
SPR->num_search_nodes = 0;
......@@ -221,17 +223,19 @@ int main(int argc, char const *argv[])
SPR->quality.function = aspectRatio;
SPR->quality.userData = NULL;
if(hxtSPR(SPR)!=0) {
int err = hxtSPR(SPR);
if(err) {
fprintf(stderr, "SPR did not work. %lu search node\n", SPR->num_search_nodes);
return 1;
}
if(SPR->tetrahedra.num==0) {
else if(SPR->tetrahedra.num==0) {
printf("there exists no tetrahedralization of that cavity\n");
return 0;
}
else {
printf("cavity tetrahedralization is outputted to out.msh\n");
return gmshTetDraw(SPR, "out.msh");
printf("cavity tetrahedralization (quality %f) is outputted to out.msh\n", SPR->tetrahedra.quality);
err = gmshTetDraw(SPR, "out.msh");
}
free(SPR);
return err;
}
\ 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