Commit 62b5fa26 authored by Célestin Marot's avatar Célestin Marot
Browse files

standalone non-recursive version... there is a bug ATM

parents
This diff is collapsed.
CFLAGS=-Wall -O3 -g -DDEBUG
LDFLAGS=-lm
EXEC=./exec
SRC= $(wildcard *.c)
OBJ= $(SRC:.c=.o)
all: $(EXEC)
$(EXEC): $(OBJ)
$(CC) -o $@ $^ $(LDFLAGS)
predicates.o: predicates.c predicates.h
hxt_SPR.o: hxt_SPR.c hxt_SPR.h predicates.h
test.o: test.c hxt_SPR.h predicates.h
$(OBJ):
# @echo $@ depends on $^
$(CC) -o $@ -c $< $(CFLAGS) -lm
run: all
@$(PWD)/$(EXEC)
clean:
rm -rf *.o
This diff is collapsed.
/* This file is part of HXTSPR. *
*
HXTSPR is free software: you can redistribute it and/or modify *
it under the terms of the GNU General Public License as published by *
the Free Software Foundation, either version 3 of the License, or *
(at your option) any later version. *
*
HXTSPR 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 General Public License for more details. *
*
You should have received a copy of the GNU General Public License *
along with HXTSPR. If not, see <http://www.gnu.org/licenses/>. *
*
See the COPYING file for the GNU General Public License . *
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#ifndef _HXT_SPR_
#define _HXT_SPR_
#ifdef __cplusplus
extern "C" {
#endif
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include "predicates.h"
// defines to set/unset the main optimizations
#define SPR_SORT_CANDIDATES
#define SPR_SAVE_QUALITIES
#define SPR_USE_FACE_HEURISTIC
// #define SPR_SAVE_CHIROTOPES
// #define SPR_INPUT_VERIFICATION
#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)
#define SPR_MAX_TETS ((SPR_MAX_PTS*SPR_MAX_PTS - 3*SPR_MAX_PTS - 2)/2)
#define SPR_MAX_FACES (3*SPR_MAX_TETS)
#if SPR_MAX_PTS > 99
#warning "Using SPR with cavities containing more" \
"than a hundred points might be prohibitively slow"
#if SPR_MAX_PTS > 255
#error "Points are indexed with uint8_t, their number is thus limited to 255."
#endif
#endif
/* structure for the points (32 bytes)*/
typedef struct {
double coord[3]; // the coordinates of the vertex (unchanged by hxtSPR)
uint32_t userID; // a user identifier for the vertex (unchanged by hxtSPR)
uint16_t valence; // used internally by the SPR function
uint8_t is_interior; // used internally by the SPR function
uint8_t is_optional; // not implemented yet (marks optional interior pts)
} SPRPoint;
/* structure for the interior, constrained edges (2 bytes)*/
typedef struct {
uint8_t node[2];
} SPREdge;
/* structure for the boundary triangles and the constrained triangles (4 bytes)
* the fourth node is used internally by the SPR function */
typedef struct {
uint8_t node[4];
} SPRTriangle;
/* a tet has the same structure as a face (4 bytes),
* except the fourth node is more useful to you */
typedef SPRTriangle SPRTet;
/* a simple bounding box structure (64 bytes) */
typedef struct {
double min[3];
double max[3];
} SPRBbox;
/* internal structure for a single step of the SPR function */
typedef struct {
uint8_t candidate[SPR_MAX_PTS-3]; // candidate points
SPRTet tet; // base face or entire tet if created
int numCandidates;
double quality; // quality before inserting this tet
} SPRStep;
/* This is the main structure defining your cavity.
* You might want to allocate this heavy structure on the heap ! */
typedef struct {
/* set max_search_nodes to limit the number of nodes in the search tree
* (as explained in the comments above the prototype of hxtSPR() )
* num_search_nodes will be set to the number of tree nodes visited by the SPR
* algorithm */
uint64_t max_search_nodes;
uint64_t num_search_nodes;
/* Points on boundary + interiors points (all points)
* hxtSPR() might shuffle what you put in here */
struct {
SPRPoint array[SPR_MAX_PTS];
int num;
} points;
/* Constrained edges
* hxtSPR() might shuffle what you put in here */
struct {
SPREdge array[SPR_MAX_EDGES];
int num;
} edges;
/* Boundary triangles + constrained triangles (appearing twice)
*
* Nodes of boundary triangles must appear in a counter-clockwise order
* when seen from the interior of the cavity.
* boundary triangles must form closed a closed surface.
* 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 */
struct {
SPRTriangle array[SPR_MAX_FACES];
SPRBbox bbox[SPR_MAX_FACES]; // filled internally
int num;
} triangles;
/* After the SPR, `tetrahedra.array` will contain the `tetrahedra.num`
* tetrahedra of the best tetrahedralization found, or be left unchanged
* if no tetrahedralization has a worse tet. with a quality above
* `tetrahedra.quality`. Note that `tetrahedra.quality` must be greater
* than or equal to zero, and it will be updated by hxtSPR to contain
* the quality of the tetrahedralization (the quality of its worse tet.)
*
* The tetrahedra will be oriented such that orient3d(a, b, c, d) < 0
* for a tet `SPRTet tet = {{a,b,c,d}};` */
struct {
SPRTet array[SPR_MAX_TETS];
int num;
double quality; // known lower bound on quality or 0
} tetrahedra;
// quality function and lower bound on quality
struct {
double (*function)(double* p0, double* p1, double* p2, double* p3,
void* userData);
void* userData;
} quality;
/* up to this point, you should not fill the structure yourself
* these are all internal stuffs */
struct {
SPRStep* array;
int num;
} steps;
struct {
#ifdef SPR_SAVE_QUALITIES
double* qualities;
#endif
#ifdef SPR_SAVE_CHIROTOPES
uint32_t* chirotope;
#endif
uint16_t* faces;
} map;
SPRExact values;
} SPRCavity;
/* hxtSPR fills the cavity with its best possible tetrahedralization in the
* sense that the quality of the worst tetrahedron is maximized.
*
* At its core, hxtSPR uses a branch and bound algorithm explained in the paper
* "Optimal Tetrahedralization for Small Polyhedron: A New Local Transformation
* 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
* 28th IMR called "Reviving the Search for Optimal Tetrahedralizations"
* by Marot Célestin, Verhetsel Kilian and Jean-François Remacle.
*
* Generally, the healthier the cavity, the faster hxtSPR will return.
* By a "healthy" cavity, we mean a cavity with a lot of possible
* tetrahedralizations and multiple high-quality ones.
*
* SPR will typically find a tetrahedralization very quickly if one exists,
* However, there may not exist any tetrahedralization, or the cavity may
* only have one very poor tetrahedralization (with lots of poor tet.)
* that the SPR will not find quickly.
*
* 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
*
* To call hxtSPR, fill the cavity structure as described in the structure
* declaration, and, if there is no error (return value is 0),
* get the resulting tetrahedralization in cavity->tetrahedra
*
*
*/
int hxtSPR(SPRCavity* cavity);
#ifdef __cplusplus
}
#endif
#endif
This diff is collapsed.
#ifndef _ROBUST_PREDICATES_H_
#define _ROBUST_PREDICATES_H_
#ifdef __cplusplus
extern "C" {
#endif
/* On some machines, the exact arithmetic routines might be defeated by the */
/* use of internal extended precision floating-point registers. Sometimes */
/* this problem can be fixed by defining certain values to be volatile, */
/* thus forcing them to be stored to memory and rounded off. This isn't */
/* a great solution, though, as it slows the arithmetic down. */
/* */
/* To try this out, write "#define INEXACT volatile" below. Normally, */
/* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */
#define INEXACT /* Nothing */
/* #define INEXACT volatile */
#define REAL double /* float or double */
/* Which of the following two methods of finding the absolute values is */
/* fastest is compiler-dependent. A few compilers can inline and optimize */
/* the fabs() call; but most will incur the overhead of a function call, */
/* which is disastrously slow. A faster way on IEEE machines might be to */
/* mask the appropriate bit, but that's difficult to do in C. */
// #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
#define Absolute(a) fabs(a)
/* we used a structure instead of the global variables used originally */
typedef struct {
/* splitter = 2^ceiling(p / 2) + 1. Used to split floats in half. */
REAL splitter;
/* A set of coefficients used to calculate maximum roundoff errors. */
REAL resulterrbound;
REAL o3derrboundB, o3derrboundC;
// Static filters for orient3d() and insphere().
// Added by H. Si, 2012-08-23.
REAL o3dstaticfilter, o3derrboundA;
} SPRExact;
void exactinit(SPRExact* values, REAL maxx, REAL maxy, REAL maxz);
REAL orient3d(const SPRExact* const values,
const REAL* const __restrict__ pa,
const REAL* const __restrict__ pb,
const REAL* const __restrict__ pc,
const REAL* const __restrict__ pd);
// remember to call exactinit with the right ranges beforehand !!!
static inline int orient3d_inline(const SPRExact* const values,
const REAL* const __restrict__ pa,
const REAL* const __restrict__ pb,
const REAL* const __restrict__ pc,
const REAL* const __restrict__ pd) {
REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
REAL det;
adx = pa[0] - pd[0];
bdx = pb[0] - pd[0];
cdx = pc[0] - pd[0];
ady = pa[1] - pd[1];
bdy = pb[1] - pd[1];
cdy = pc[1] - pd[1];
adz = pa[2] - pd[2];
bdz = pb[2] - pd[2];
cdz = pc[2] - pd[2];
bdxcdy = bdx * cdy;
cdxbdy = cdx * bdy;
cdxady = cdx * ady;
adxcdy = adx * cdy;
adxbdy = adx * bdy;
bdxady = bdx * ady;
det = adz * (bdxcdy - cdxbdy)
+ bdz * (cdxady - adxcdy)
+ cdz * (adxbdy - bdxady);
int ret = (det > values->o3dstaticfilter) - (det < -values->o3dstaticfilter);
if (ret!=0) return ret;
det = orient3d(values, pa,pb,pc,pd);
return (det>0.0) - (det<0.0);
}
#ifdef __cplusplus
}
#endif
#endif
/* This file is part of HXTSPR. *
*
HXTSPR is free software: you can redistribute it and/or modify *
it under the terms of the GNU General Public License as published by *
the Free Software Foundation, either version 3 of the License, or *
(at your option) any later version. *
*
HXTSPR 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 General Public License for more details. *
*
You should have received a copy of the GNU General Public License *
along with HXTSPR. If not, see <http://www.gnu.org/licenses/>. *
*
See the COPYING file for the GNU General Public License . *
*
Author: Célestin Marot (celestin.marot@uclouvain.be) */
#include "hxt_SPR.h"
#include <math.h>
double aspectRatio(double a[3], double b[3], double c[3], double d[3],
void* userData)
{
(void)(userData); // unused
double ab[3], ac[3], ad[3], bc[3], cd[3], db[3];
for (int i=0; i<3; i++) {
ab[i] = b[i] - a[i]; // AB
ac[i] = c[i] - a[i]; // AC
ad[i] = d[i] - a[i]; // AD
}
double acxad0 = ac[1]*ad[2] - ac[2]*ad[1];
double adxab0 = ad[1]*ab[2] - ad[2]*ab[1];
double abxac0 = ab[1]*ac[2] - ab[2]*ac[1];
double volume6 = ab[0]*acxad0 + ac[0]*adxab0 + ad[0]*abxac0;
// abort as early as possible
if(volume6<=0.0)
return 0.0;
double acxad1 = ac[2]*ad[0] - ac[0]*ad[2];
double acxad2 = ac[0]*ad[1] - ac[1]*ad[0];
double adxab1 = ad[2]*ab[0] - ad[0]*ab[2];
double adxab2 = ad[0]*ab[1] - ad[1]*ab[0];
double abxac1 = ab[2]*ac[0] - ab[0]*ac[2];
double abxac2 = ab[0]*ac[1] - ab[1]*ac[0];
for (int i=0; i<3; i++) {
db[i] = b[i] - d[i]; // DB = B-D = AB-AD
bc[i] = c[i] - b[i]; // BC = C-B = AC-AB
cd[i] = d[i] - c[i]; // CD = D-c = AD-AC
}
double bcxcd0 = bc[1]*cd[2] - bc[2]*cd[1]; // = acxad0+abxac0+adxab0;
double bcxcd1 = bc[2]*cd[0] - bc[0]*cd[2]; // = acxad1+abxac1+adxab1;
double bcxcd2 = bc[0]*cd[1] - bc[1]*cd[0]; // = acxad2+abxac2+adxab2;
double areaSum = sqrt(acxad0*acxad0 + acxad1*acxad1 + acxad2*acxad2)
+ sqrt(adxab0*adxab0 + adxab1*adxab1 + adxab2*adxab2)
+ sqrt(abxac0*abxac0 + abxac1*abxac1 + abxac2*abxac2)
+ sqrt(bcxcd0*bcxcd0 + bcxcd1*bcxcd1 + bcxcd2*bcxcd2);
double l = ab[0]*ab[0] + ab[1]*ab[1] + ab[2]*ab[2]; // |AB|²
double l2 = ac[0]*ac[0] + ac[1]*ac[1] + ac[2]*ac[2]; // |AC|²
double l3 = ad[0]*ad[0] + ad[1]*ad[1] + ad[2]*ad[2]; // |AD|²
double l4 = bc[0]*bc[0] + bc[1]*bc[1] + bc[2]*bc[2]; // |BC|²
double l5 = cd[0]*cd[0] + cd[1]*cd[1] + cd[2]*cd[2]; // |CD|²
double l6 = db[0]*db[0] + db[1]*db[1] + db[2]*db[2]; // |DB|²
if(l2>l) l=l2;
if(l3>l) l=l3;
if(l4>l) l=l4;
if(l5>l) l=l5;
if(l6>l) l=l6;
return sqrt(24*volume6*volume6/(l*areaSum*areaSum));
}
// super simple visualisation with gmsh
int gmshTetDraw(SPRCavity* SPR, const char* filename)
{
FILE* file = fopen(filename,"w");
if(file==NULL)
return 1;
/* format for gmsh */
fprintf(file,"$MeshFormat\n"
"2.2 0 %u\n"
"$EndMeshFormat\n"
"$Nodes\n"
"%d\n",(unsigned) sizeof(double), SPR->points.num);
for (int i=0; i<SPR->points.num; i++)
fprintf(file,"%u %.10E %.10E %.10E\n",i+1, SPR->points.array[i].coord[0],
SPR->points.array[i].coord[1],
SPR->points.array[i].coord[2]);
fprintf(file,"$EndNodes\n"
"$Elements\n"
"%d\n", SPR->tetrahedra.num);
for (int i=0; i<SPR->tetrahedra.num; i++){
fprintf(file,"%d 4 0 %u %u %u %u\n", i+1, SPR->tetrahedra.array[i].node[0]+1,
SPR->tetrahedra.array[i].node[1]+1,
SPR->tetrahedra.array[i].node[2]+1,
SPR->tetrahedra.array[i].node[3]+1);
}
fputs("$EndElements\n",file);
fclose(file);
return 0;
}
int main(int argc, char const *argv[])
{
double coord[25][3] = {
{0.58671240234164312443709832223248668015003204345703,0.03516283877473333485674800158449215814471244812012,0.04997715958560345478334241420270700473338365554810},
{0.69999999999999995559107901499373838305473327636719,-0.00000000000000017145055188062938983975907254829012,-0.00000000000000004898587196589413075214088999399776},
{0.67320508075688778593104188985307700932025909423828,-0.00000000000000016488768946373159238004406226972037,-0.09999999999999993616217608405349892564117908477783},
{0.57955393645644781575043680277303792536258697509766,0.07245318826589322580566943088342668488621711730957,-0.05343877417631059945080096440506167709827423095703},
{0.68234535975465848700594051479129120707511901855469,0.07450181889388610145807945173146435990929603576660,-0.07249024801756982772094062283940729685127735137939},
{0.53946158431996149573706134106032550334930419921875,-0.01967558269480869370671172191578079946339130401611,-0.08238916464139039996794622311426792293787002563477},
{0.69179629664624397999972416073433123528957366943359,0.10685449895321200153297525048401439562439918518066,-0.00000000000000004898587196589413075214088999399776},
{0.56780156011846882879723352743894793093204498291016,-0.06076576841659741406997952140045526903122663497925,0.00923681616527150493867193148389560519717633724213},
{0.63797938422767774824251318932510912418365478515625,0.08734851713085420299442773739428957924246788024902,-0.13886610280937750072638436904526315629482269287109},
{0.49184110864420532704954780456318985670804977416992,0.09270955818833073736051630930887768045067787170410,-0.00024693681951923340811561047303257510066032409668},
{0.60000000000000019984014443252817727625370025634766,-0.00000000000000014695761589768240458237431406030273,-0.17320508075688759164201258045068243518471717834473},
{0.68201852245367466931469380142516456544399261474609,0.07516371438559803219714439137533190660178661346436,0.07313675010055335767855666517789359204471111297607},
{0.56157280600954107008249138743849471211433410644531,0.13340233585380706360368208152067381888628005981445,0.08058263496460443819913166407786775380373001098633},
{0.68657610052058659633189563464839011430740356445312,-0.06580504799977010543798883190902415663003921508789,-0.06328819191132203925143073774961521849036216735840},
{0.68192597680842104157505900730029679834842681884766,-0.07447818796802643803722787652077386155724525451660,0.07355976420064183618485742499615298584103584289551},
{0.44625545218884238307666123546368908137083053588867,-0.02013109738830444819801535061287722783163189888000,-0.03277409574633017436129378552323032636195421218872},
{0.61506778805773154594760399049846455454826354980469,-0.07370695677550258284593809321449953131377696990967,-0.16039732307644941067259480860229814425110816955566},
{0.67320508075688767490873942733742296695709228515625,-0.00000000000000016488768946373159238004406226972037,0.10000000000000019984014443252817727625370025634766},
{0.69179629664624386897742169821867719292640686035156,-0.10685449895321259827785098650565487332642078399658,-0.00000000000000004898587196589413075214088999399776},
{0.63813593530542289666840360951027832925319671630859,-0.08788260128750351174442556612120824865996837615967,0.13862965526503509550160231356130680069327354431152},
{0.49173216365874278244163519957510288804769515991211,-0.01984501837771410898136892342336068395525217056274,0.09065199313269833858441870688693597912788391113281},
{0.65601444577662237911397369316546246409416198730469,-0.11549518170588020382627547633092035539448261260986,-0.11139826339487800266336847698767087422311305999756},
{0.63825557951601918915685018873773515224456787109375,0.08962706452937579948425650400167796760797500610352,0.13825557951601941120145511376904323697090148925781},
{0.54864718139337009805700517972582019865512847900391,0.11577765142005609633191198781787534244358539581299,0.19055671192273110348125442214950453490018844604492},
{0.59999999999999975575093458246556110680103302001953,-0.00000000000000014695761589768230597476116143382705,0.17320508075688789695334435236873105168342590332031}
};
unsigned char face[44][3] = {
{2, 1, 4},
{4, 1, 6},
{3, 4, 6},
{8, 2, 4},
{3, 8, 4},
{9, 5, 3},
{3, 5, 10},
{6, 1, 11},
{12, 3, 6},
{1, 2, 13},
{10, 2, 8},
{3, 10, 8},
{15, 5, 9},
{12, 9, 3},
{2, 10, 16},
{16, 10, 5},
{15, 7, 5},
{5, 7, 13},
{11, 1, 17},
{12, 6, 11},
{13, 18, 1},
{18, 13, 7},
{1, 18, 14},
{14, 18, 7},
{14, 17, 1},
{14, 7, 19},
{9, 20, 15},
{20, 9, 12},
{16, 5, 21},
{21, 2, 16},
{20, 7, 15},
{5, 13, 21},
{13, 2, 21},
{17, 22, 11},
{11, 22, 12},
{14, 19, 17},
{19, 7, 20},
{23, 20, 12},
{17, 24, 22},
{22, 23, 12},
{24, 17, 19},
{20, 24, 19},
{23, 24, 20},
{22, 24, 23}
};
SPRCavity* SPR = malloc(sizeof(SPRCavity));
SPR->max_search_nodes = UINT64_MAX;
SPR->num_search_nodes = 0;
SPR->points.num = sizeof(coord)/sizeof(coord[0]);
for (int i=0; i<SPR->points.num; i++) {
for (int j=0; j<3; j++)
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->edges.num = 0;
SPR->tetrahedra.num = 0;
SPR->tetrahedra.quality = 0.0;
// SPR->tetrahedra.quality = 0.27985993678969944; // giving a lower bound for the quality
SPR->quality.function = aspectRatio;
SPR->quality.userData = NULL;
if(hxtSPR(SPR)!=0) {
fprintf(stderr, "SPR did not work. %lu search node\n", SPR->num_search_nodes);
return 1;
}
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");
}
}
\ 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