summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorUrbain Vaes <u.vaes13@imperial.ac.uk>2015-04-25 18:40:00 +0100
committerUrbain Vaes <u.vaes13@imperial.ac.uk>2015-04-25 18:40:00 +0100
commit1d624269facbde200d214e27a1f6d678fbd01ab5 (patch)
tree2a16e04755150c7590ef9454e2095ebcf59654b5
parent7fbbe36d7da913729660d0987c09050f152248cf (diff)
Upload of the project
-rw-r--r--ExteriorNodes.c106
-rw-r--r--Makefile6
-rw-r--r--ShowIt.c234
-rw-r--r--Split.c220
-rw-r--r--addNode.c540
-rw-r--r--addTriangle.c75
-rw-r--r--checkAndSplitEncroachedSegments.c281
-rw-r--r--checkAndSplitSkinnyTriangle.c164
-rw-r--r--freeStructures.c67
-rw-r--r--header.h130
-rw-r--r--initStructures.c106
-rw-r--r--isInFunctions.c306
-rw-r--r--main.c175
-rw-r--r--makefile.linux7
-rw-r--r--printManPage.c22
-rw-r--r--readInputFile.c150
-rw-r--r--removeExteriorTriangles.c39
-rw-r--r--superior.poly1045
18 files changed, 3673 insertions, 0 deletions
diff --git a/ExteriorNodes.c b/ExteriorNodes.c
new file mode 100644
index 0000000..46799f0
--- /dev/null
+++ b/ExteriorNodes.c
@@ -0,0 +1,106 @@
+#include "header.h"
+
+/*
+OBJECTIFS:
+addExteriorNodes est une fonction permettant d'ajouter 4 points englobant le
+nuage de points lu par le programme. Le but est de faciliter l'ajout des points
+dans la triangulation en s'assurant que chaque nouveau point sera situe a
+l'interieur d'un triangle existant. Les 4 points ajoutes doivent etre suffisamment
+eloignes des autres points afin de conserver la convexite de la forme a triangulariser.
+
+INPUTS:
+-> Structure theGrid contenant les coordonnees de l'ensemble des points fournis.
+-> Deux vecteurs seg1 et seg2 qui contiennent la liste des extremites des segments.
+-> Structure theTriangulation regroupant les pointeurs vers les triangles de la
+ triangulation courante (initialement vide).
+
+OUTPUTS:
+-> La structure theGrid est mise à jour et contient les 4 nouveaux points en bout
+ de liste.
+-> La triangulation est remplie avec deux premiers triangles, traces a partir des
+ 4 points exterieurs.
+-> Les triangles sont ajoutés dans l'arbre qui trie les triangles suivant leur qualite.
+-> 4 nouveaux segments (les 4 cotes du grand carre) sont ajoutes a la liste des segments.
+*/
+
+void addExteriorNodes(structGrid *theGrid,structVector *seg1,structVector *seg2,structTriangulation *theTriangulation)
+{
+ // Etape 1 : recherche des points extremes du nuages de points fournis
+
+ double minX = theGrid->X[0];
+ double maxX = theGrid->X[0];
+ double minY = theGrid->Y[0];
+ double maxY = theGrid->Y[0];
+
+ int n = theGrid->nNode;
+ theGrid->nNode = theGrid->nNode+4;
+ theGrid->nInputVertex = theGrid->nNode;
+
+ int i;
+
+ for(i=1;i<n;i++)
+ {
+ minX = Min(minX,theGrid->X[i]);
+ maxX = Max(maxX,theGrid->X[i]);
+ minY = Min(minY,theGrid->Y[i]);
+ maxY = Max(maxY,theGrid->Y[i]);
+ }
+
+ // Etape 2 : definition des 4 points exterieurs. Le parametre "aa" permet
+ // d'eloigner plus ou moins les 4 points exterieurs des autres points
+
+ double aa = 0.2;
+ theGrid->X[n] = minX-(maxX-minX)*aa;
+ theGrid->Y[n] = minY-(maxY-minY)*aa;
+
+ theGrid->X[n+1] = maxX+(maxX-minX)*aa;
+ theGrid->Y[n+1] = minY-(maxY-minY)*aa;
+
+ theGrid->X[n+2] = minX-(maxX-minX)*aa;
+ theGrid->Y[n+2] = maxY+(maxY-minY)*aa;
+
+ theGrid->X[n+3] = maxX+(maxX-minX)*aa;
+ theGrid->Y[n+3] = maxY+(maxY-minY)*aa;
+
+ // Etape 3 : defintion de deux triangles et ajout dans la triangulation
+
+ structTriangle *Triangle1 = TriangleNew(n+2,n+1,n+3,theGrid);
+ structTriangle *Triangle2 = TriangleNew(n+2,n+3,n+4,theGrid);
+
+ Triangle1->index = 0;
+ Triangle2->index = 1;
+
+ Triangle1->Voisins[0] = -1;
+ Triangle1->Voisins[1] = -2;
+ Triangle1->Voisins[2] = Triangle2->index;
+
+ Triangle2->Voisins[0] = Triangle1->index;
+ Triangle2->Voisins[1] = -3;
+ Triangle2->Voisins[2] = -4;
+
+ theTriangulation->Triangles[Triangle1->index] = Triangle1;
+ theTriangulation->Triangles[Triangle2->index] = Triangle2;
+ theTriangulation->nTriangles = 2;
+
+
+ // Etape 4 : creation des 4 nouveaux segments
+
+ int m = seg1->nElem;
+ seg1->nElem = seg1->nElem + 4;
+ seg2->nElem = seg2->nElem + 4;
+
+ seg1->elem[m] = n+1; seg2->elem[m] = n+2;
+ seg1->elem[m+1] = n+2; seg2->elem[m+1] = n+4;
+ seg1->elem[m+2] = n+4; seg2->elem[m+2] = n+3;
+ seg1->elem[m+3] = n+3; seg2->elem[m+3] = n+1;
+
+ // Etape 5 : les triangles sont ajoutes dans l'arbre de tri
+
+ theTriangulation->tree->triangle = theTriangulation->Triangles[0];
+ theTriangulation->tree->quality = theTriangulation->Triangles[0]->quality;
+ theTriangulation->tree->left = QualityTreeNew();
+ theTriangulation->tree->right = QualityTreeNew();
+ theTriangulation->begin = theTriangulation->tree;
+
+}
+
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..00e44f7
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,6 @@
+COMPILER = gcc -O3
+OBJECTS = *.c
+OUTPUT = ruppert
+
+$(OUTPUT) : $(OBJECTS)
+ $(COMPILER) -framework Cocoa -framework OpenGL -lglfw -lm $(OBJECTS) -o $(OUTPUT)
diff --git a/ShowIt.c b/ShowIt.c
new file mode 100644
index 0000000..f287135
--- /dev/null
+++ b/ShowIt.c
@@ -0,0 +1,234 @@
+#include "header.h"
+
+const float rotations_per_tick = .2;
+const float zoom_per_tick = 1;
+const float translate_per_tick = 1;
+
+void showIt(structGrid *theGrid, structTriangulation *theTriangulation,structVector *seg1,structVector *seg2, int Dimension2, int withCircles,double angle)
+{
+ float rotate_y = 180;
+ float rotate_z = 0;
+ float rotate_x;
+ if(Dimension2==0){rotate_x = -60.0;}
+ else {rotate_x = 0.0;}
+ float translate_x = 0;
+ float translate_y = 0;
+
+ double* x = theGrid->X;
+ double* y = theGrid->Y;
+
+ const int window_width = 800,
+ window_height = 600;
+
+ if (glfwInit() != GL_TRUE)
+ Shut_Down(1);
+
+ float translate_z = ReshapeAndShift(theGrid);
+
+ if (glfwOpenWindow(window_width, window_height, 5, 6, 5,
+ 0, 0, 0, GLFW_WINDOW) != GL_TRUE)
+ Shut_Down(1);
+ glfwSetWindowTitle("Geometrie numerique en CAO - projet (partie2) - Groupe 6");
+
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+ float aspect_ratio = ((float)window_height) / window_width;
+ double a = .5;
+ glFrustum(a, -a, -a * aspect_ratio, a * aspect_ratio, 1, 100000);
+ glMatrixMode(GL_MODELVIEW);
+
+
+ double old_time = glfwGetTime();
+ while(1)
+ {
+ double current_time = glfwGetTime(),
+ delta_rotate = (current_time - old_time) * rotations_per_tick * 360,
+ delta_zoom = (current_time - old_time) * zoom_per_tick,
+ delta_translate = (current_time - old_time) *translate_per_tick*translate_z;
+
+ old_time = current_time;
+
+ if (glfwGetKey(GLFW_KEY_ESC) == GLFW_PRESS)
+ break;
+ if (glfwGetKey(GLFW_KEY_LEFT) == GLFW_PRESS)
+ rotate_z += delta_rotate;
+ if (glfwGetKey(GLFW_KEY_RIGHT) == GLFW_PRESS)
+ rotate_z -= delta_rotate;
+
+ if (glfwGetKey(GLFW_KEY_UP) == GLFW_PRESS)
+ rotate_x -= delta_rotate;
+ if (glfwGetKey(GLFW_KEY_DOWN) == GLFW_PRESS)
+ rotate_x += delta_rotate;
+
+ if (glfwGetKey('T') == GLFW_PRESS)
+ translate_z *= (1-delta_zoom);
+
+ if (glfwGetKey('R') == GLFW_PRESS)
+ translate_z *= 1+delta_zoom;
+
+ if (glfwGetKey('O') == GLFW_PRESS)
+ translate_y += delta_translate;
+ if (glfwGetKey('L') == GLFW_PRESS)
+ translate_y -= delta_translate;
+
+ if (glfwGetKey('I') == GLFW_PRESS)
+ translate_x += delta_translate;
+ if (glfwGetKey('P') == GLFW_PRESS)
+ translate_x -= delta_translate;
+
+ glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+ glLoadIdentity();
+ glTranslatef(translate_x, translate_y, translate_z);
+
+ glRotatef(rotate_x, 1, 0, 0);
+ glRotatef(rotate_z, 0, 0, 1);
+ glRotatef(rotate_y, 0, 1, 0);
+
+ double* z = theGrid->Z;
+ int nN = theGrid->nNode;
+ int nT = theTriangulation->nTriangles;
+ Draw_Triangulation(nN,x,y,z,theTriangulation,nT,Dimension2,withCircles,seg1,seg2,angle);
+
+ glfwSwapBuffers();
+
+ }
+ Shut_Down(0);
+}
+
+
+
+void Shut_Down(int return_code)
+{
+ glfwTerminate();
+ //exit(return_code);
+}
+
+void Draw_Triangulation(int nNodes, double *x, double *y, double *z, structTriangulation *T, int nTriang,int Dimension2, int withCircles,structVector *seg1,structVector *seg2,double angle)
+{
+ // Etape 0 : defition de la qualite des triangles
+
+ double minAngle = angle*3.1415/180.0;
+ double qualityRequest = 1.0/(2.0*sin(minAngle));
+
+ // Etape 1 : Dans le cas 2D, les composantes "z" sont mises a zero
+
+ int i; double aa;
+ if(Dimension2==1)
+ {
+ aa = 0.0;
+ }
+ else
+ {
+ aa = 1.0;
+ }
+
+ // Etape 2 : Dessin des aretes des triangles. Les segments sont dessines
+ // en vert, les triangles de mauvaise qualite sont dessines en rouge et
+ // les triangles de bonne qualite sont dessines en bleu.
+
+ for(i=0;i<nTriang;++i)
+ {
+ structTriangle *Triang = T->Triangles[i];
+
+ if(Triang->toBeDeleted==0)
+ {
+ glLineWidth(1.0);
+ glPointSize(4.0);
+
+ glBegin(GL_LINES);
+ {
+ int j;
+ for(j=0;j<3;++j)
+ {
+ int V1 = Triang->Nodes[j];
+ int V2 = Triang->Nodes[(j+1)%3];
+ if(Triang->Segment[j] != -1)
+ glColor3f(0, 1, 0);
+ else
+ glColor3f(0, 0, 1);
+
+ if(Triang->quality > 1)
+ glColor3f(1, 0, 0);
+
+ glVertex3f(x[V1-1],y[V1-1],aa*z[V1-1]);
+ glVertex3f(x[V2-1],y[V2-1],aa*z[V2-1]);
+
+ }
+ }
+ glEnd();
+ }
+ }
+
+ // Etape 3 : Si l'utilisateur l'a demande, les cercles circonscrits aux triangles
+ // seront dessines.
+
+ if(withCircles==1)
+ {
+ for(i=0;i<nTriang;++i)
+ {
+ structTriangle *Triang = T->Triangles[i];
+ if(Triang->toBeDeleted==0)
+ {
+ glPointSize(1.0);
+ DrawCircle(Triang, x, y);
+ }
+ }
+ }
+}
+
+
+void DrawCircle(structTriangle *Triangle, double *X, double *Y)
+{
+ structGrid *theGrid = malloc(sizeof(structGrid));
+ theGrid->X = X;
+ theGrid->Y = Y;
+ double CenterX,CenterY;
+ computeCircumcenter(Triangle,theGrid,&CenterX,&CenterY);
+ double x = X[Triangle->Nodes[0]-1];
+ double y = Y[Triangle->Nodes[0]-1];
+ double radius = sqrt((x-CenterX)*(x-CenterX)+(y-CenterY)*(y-CenterY));
+
+ glColor3f (1.0, 0.0, 0.0);
+ const float TWOPI = 6.283185;
+ int slices = 50;
+ float step = TWOPI/slices;
+ glBegin(GL_LINE_LOOP);
+ float angle = 0.0;
+ for(angle = 0.0; angle < TWOPI; angle+=step)
+ glVertex2f(CenterX+cos(angle) * radius, CenterY + sin(angle) * radius);
+ glEnd();
+}
+
+float ReshapeAndShift(structGrid *theGrid)
+{
+ double* x = theGrid->X;
+ double* y = theGrid->Y;
+
+ double minX = min(x,theGrid->nNode);
+ double maxX = max(x,theGrid->nNode);
+ double minY = min(y,theGrid->nNode);
+ double maxY = max(y,theGrid->nNode);
+ double sizeX = (maxX-minX)/1.9;
+ double sizeY = (maxY-minY)/1.9;
+
+ double size = Max(sizeX,sizeY);
+
+ return -size*3;
+}
+
+double min(double *x, int n) {
+ double myMin = x[0];
+ int i;
+ for (i=1 ;i < n; i++)
+ if (x[i] < myMin) myMin = x[i];
+ return myMin;
+}
+
+double max(double *x, int n) {
+ double myMax = x[0];
+ int i;
+ for (i=1 ;i < n; i++)
+ if (x[i] > myMax) myMax = x[i];
+ return myMax;
+}
diff --git a/Split.c b/Split.c
new file mode 100644
index 0000000..d94f39d
--- /dev/null
+++ b/Split.c
@@ -0,0 +1,220 @@
+#include "header.h"
+
+/*
+OBJECTIFS:
+Lorsqu'un segment est empiete par un noeud, cette fonction calcule les coordonnees du
+point de division. Plusieurs cas sont a considerer:
+-> Lorsque le noeud qui empiete sur le segment n'appartient pas a un autre segment
+ alors le segment est divise en son milieu.
+-> Lorsque aucune extremite du segment a diviser ne fait partie des noeuds fournis
+ dans le fichier d'input, alors le segement est divise en son milieu (case 0).
+-> Lorsque les deux extremites du segment a diviser font partie des noeuds fournis
+ dans le fichier d'input, alors le segement est divise en son milieu (case 2).
+-> Dans les autres cas, le segment est divise selon le principe des "concentric shells":
+ la distance entre l'extremite du segment et le point de division est une puissance
+ de 2 d'une longueur D fixee a 0.1% de la longueur du plus petit segment (case 1).
+
+
+INPUTS:
+-> theGrid : coordonnees de l'ensemble des points du maillage.
+-> seg1 et seg2 : liste des extremites des segments.
+-> i: numero du segment a diviser.
+-> SegPointEncroaching: vaut 1 si le segment est empiete par un point qui se trouve sur
+ un autre segment.
+
+OUTPUTS:
+-> x et y contiennent les coordonnees du point auquel sera divise le segment i.
+*/
+void setSplittingPoint(structGrid *theGrid, structVector *seg1, structVector *seg2, int i, double *x, double *y, int SegPointEncroaching)
+{
+ // Etape 1 : coordonnees des extremites du segment
+
+ double ax = theGrid->X[seg1->elem[i]-1];
+ double ay = theGrid->Y[seg1->elem[i]-1];
+ double bx = theGrid->X[seg2->elem[i]-1];
+ double by = theGrid->Y[seg2->elem[i]-1];
+
+
+ // Etape 2 : recherche de la configuration
+
+ int theCase=0, a=0,b=0;
+
+ if(seg1->elem[i] <= theGrid->nInputVertex)
+ {
+ theCase++;
+ a = 1;
+ }
+ if(seg2->elem[i] <= theGrid->nInputVertex)
+ {
+ theCase++;
+ b = 1;
+ }
+
+ // Etape 3 : Calcul des coordonnees du point de division du segment
+
+ double D = 0.0001*seg1->lengthMinSeg;
+
+ // Etape 3a : Division au milieu du segment
+ if(theCase == 0 || theCase == 2 || !SegPointEncroaching)
+ {
+ *x = 0.5*(ax + bx);
+ *y = 0.5*(ay + by);
+ }
+ // Etape 3b : "Concentric shells"
+ else if(theCase == 1)
+ {
+ int k,z; k = 0;
+ double d = 0.5*sqrt((ax-bx)*(ax-bx)+(ay-by)*(ay-by));
+ double dist = D;
+ while(dist<d)
+ {
+ k++;
+ dist = dist*2;
+ }
+
+ if (fabs(dist-d) < fabs(0.5*dist-d))
+ {
+ z = k;
+ dist = dist;
+ }
+ else
+ {
+ z = k-1;
+ dist = dist*0.5;
+ }
+
+ double w1 = a*(1-dist/(2*d))+b*dist/(2*d);
+ double w2 = 1-w1;
+
+ *x = w1*ax+w2*bx;
+ *y = w1*ay+w2*by;
+ }
+}
+
+
+/*
+OBJECTIFS:
+Cette fonction permet de diviser un segment en deux sous segments, lorsque ce segment
+est empiete par un point du maillage.
+
+INPUTS:
+-> T : pointeur vers un triangle adjacent au segment.
+-> theTriangulation : ensemble des triangles de la triangulation courante.
+-> theGrid : coordonnees de l'ensemble des points du maillage.
+-> seg1 et seg2 : liste des extremites des segments.
+-> i: numero du segment a diviser.
+-> SegPointEncroaching: vaut 1 si le segment est empiete par un point qui se trouve sur
+ un autre segment.
+
+OUTPUTS:
+-> Le segment a ete divise en deux parties. Le point de division a ete ajoute a la liste
+ de noeuds, la triangulation a ete mise a jour et les deux sous segments ont ete ajoutes
+ dans la liste des segments.
+*/
+void splitSeg(structTriangle *T, structTriangulation *theTriangulation, structGrid *theGrid, structVector *seg1, structVector *seg2, int i, int SegPointEncroaching)
+{
+ // Etape 1 : nouvelle allocation de memoire pour theGrid si necessaire.
+
+ int nn = theGrid->nNode;
+ if(nn>=theGrid->size)
+ {
+ int N = 2;
+ theGrid->X = realloc(theGrid->X,N*theGrid->size*sizeof(double));
+ theGrid->Y = realloc(theGrid->Y,N*theGrid->size*sizeof(double));
+ theGrid->Z = realloc(theGrid->Z,N*theGrid->size*sizeof(double));
+ theGrid->size = N*theGrid->size;
+ //printf("\n\n----------------------------------------------------- \n");
+ //printf("Realloc of theGrid->X, theGrid->Y and theGrid->Z \n");
+ //printf("new size %d - nNode %d - nInputVertex %d \n",theGrid->size,theGrid->nNode,theGrid->nInputVertex);
+ //printf("----------------------------------------------------- \n");
+ }
+
+
+ // Etape 2 : calcul des coordonnees du point de division du segement.
+
+ double x,y;
+ setSplittingPoint(theGrid, seg1, seg2, i, &x, &y, SegPointEncroaching);
+
+
+ // Etape 3 : ajout du nouveau point dans la liste de noeuds.
+
+ theGrid->X[nn] = x;
+ theGrid->Y[nn] = y;
+ theGrid->nNode++;
+
+
+ // Etape 4 : ajout des deux sous segments dans la liste des segments
+ // et nouvelle allocation de memoire si necessaire.
+
+ seg2->elem[seg2->nElem] = seg2->elem[i];
+ seg2->elem[i] = nn+1;
+ seg1->elem[seg1->nElem] = nn+1;
+
+ seg1->nElem++;
+ seg2->nElem++;
+
+ if(seg1->nElem>=seg1->size)
+ {
+ //printf("\n\n----------------------------------------------------- \n");
+ //printf("Realloc of seg1->elem and seg2->elem \n");
+ //printf("----------------------------------------------------- \n");
+ int N = 2;
+ seg1->elem = realloc(seg1->elem,N*seg1->size*sizeof(int));
+ seg2->elem = realloc(seg2->elem,N*seg2->size*sizeof(int));
+ seg1->size = seg1->size*N;
+ seg2->size = seg1->size*N;
+ }
+
+
+ // Etape 5 : mise a jour de la triangulation.
+
+ addNode(T, nn+1,theTriangulation, theGrid,seg1,seg2);
+}
+
+/*
+OBJECTIFS:
+Cette fonction permet de diviser un triangle de mauvaise qualite en trois sous triangles.
+Un nouveau noeud correspondant au centre du cercle circonscrit au triangle est ajoute.
+
+INPUTS:
+-> theTriangle : pointeur vers un triangle de mauvaise qualite.
+-> theTriangulation : ensemble des triangles de la triangulation courante.
+-> theGrid : coordonnees de l'ensemble des points du maillage.
+-> px et py : coordonnees du centre du cercle circonscrit au triangle.
+-> seg1 et seg2 : liste des extremites des segments.
+-> i: numero du segment a diviser.
+
+OUTPUTS:
+-> Le triangle a ete divise en trois parties. Le point de division a ete ajoute a la liste
+ de noeuds et la triangulation a ete mise a jour.
+*/
+void splitTri(structTriangle *theTriangle, structTriangulation *theTriangulation, structGrid *theGrid, double px, double py, structVector *seg1, structVector *seg2)
+{
+ // Etape 1 : nouvelle allocation de memoire pour theGrid si necessaire.
+
+ int nn = theGrid->nNode;
+ if(nn>=theGrid->size)
+ {
+ int N = 2;
+ theGrid->X = realloc(theGrid->X,N*theGrid->size*sizeof(double));
+ theGrid->Y = realloc(theGrid->Y,N*theGrid->size*sizeof(double));
+ theGrid->Z = realloc(theGrid->Z,N*theGrid->size*sizeof(double));
+ theGrid->size = N*theGrid->size;
+ //printf("\n\n----------------------------------------------------- \n");
+ //printf("Realloc of theGrid->X, theGrid->Y and theGrid->Z \n");
+ //printf("new size %d - nNode %d - nInputVertex %d \n",theGrid->size,theGrid->nNode,theGrid->nInputVertex);
+ //printf("----------------------------------------------------- \n");
+ }
+
+
+ // Etape 2 : ajout du nouveau point a la liste des noeuds.
+
+ theGrid->X[nn] = px;
+ theGrid->Y[nn] = py;
+ theGrid->nNode++;
+
+
+ // Etape 3 : mise a jour de la triangulation.
+ addNode(theTriangle, nn+1,theTriangulation, theGrid,seg1,seg2);
+}
+
diff --git a/addNode.c b/addNode.c
new file mode 100644
index 0000000..da94987
--- /dev/null
+++ b/addNode.c
@@ -0,0 +1,540 @@
+
+#include "header.h"
+
+/*
+OBJECTIFS:
+Permet de mettre a jour la triangulation de Delaunay suite a l'ajout d'un noeud:
+1) Determiner dans quel triangle se trouve le point d'indice indexNode (par la methode de la secante).
+2) Trouver la cavite associee a ce point.
+3) Repertorier les aretes formant le contour de cette cavite.
+4) Remplir la cavite de nouveaux triangles, dont un des sommets est indexNode.
+
+INPUTS:
+-> theTriangle : pointeur vers un triangle susceptible d'etre proche du point a ajouter.
+-> indexNode : indice du point a ajouter dans la triangulation.
+-> theMesh : strucuture contenant les pointeurs vers les triangles de la triangulation courante.
+-> theGrid: structure contenant les coordonnees des noeuds du maillage.
+-> seg1 et seg2 : liste des extremites des segments.
+
+OUTPUT:
+-> La triangulation a ete mise a jour, de nouveaux triangles ont ete crees pour contenir le
+ nouveau point. Ces nouveaux triangles ont ete ajoutes dans l'arbre qui trie les triangles
+ suivant leur qualite.
+*/
+
+void addNode(structTriangle *theTriangle, int indexNode, structTriangulation *theMesh, structGrid *theGrid, structVector *seg1, structVector *seg2)
+{
+ // Etape 1 : Nouvelle allocation de memoire pour la triangulation si necessaire.
+
+ if(theMesh->nTriangles>=theMesh->size-20)
+ {
+ //printf("\n\n----------------------------------------------------- \n");
+ //printf("Realloc of theMesh->Triangles\n");
+ //printf("----------------------------------------------------- \n");
+ int N = 2;
+ theMesh->Triangles = realloc(theMesh->Triangles,N*theMesh->size*sizeof(structTriangle));
+ theMesh->size = theMesh->size*N;
+ }
+ if(theGrid->nNode%25000 == 0){printf("Nodes: %d \n", theGrid->nNode);}
+
+
+ // Etape 2 : Recherche d'un triangle dont le cercle circonscrit contient le point
+ // a ajouter dans la triangulation.
+
+ double x = theGrid->X[indexNode-1];
+ double y = theGrid->Y[indexNode-1];
+ theTriangle = FindFirstTri(theTriangle, x, y, theGrid, theMesh,-1,1);
+
+
+ // Etape 3 : Determiner les triangles de la cavite.
+
+ structVector *Cavity = newVector(1000);
+ buildCavity(Cavity,theMesh,theTriangle,x,y,theGrid);
+
+
+ // Etape 4 : Determiner les aretes du contour de la cavite.
+
+ structVector *Edge1 = newVector(3*Cavity->nElem);
+ structVector *Edge2 = newVector(3*Cavity->nElem);
+ GetEdges(Edge1, Edge2, Cavity, theMesh);
+
+
+ // Etape 5 : Recuperer les indices des triangles de la cavite afin d'utiliser ulterieurement
+ // ces nouveaux indices pour stocker les nouveaux triangles dans la tringulation. La surface
+ // de la cavite est egalement calculee.
+
+ int k; int n = Edge1->nElem;
+ int* newIndexes = malloc(sizeof(int)*n);
+ double AreaCavity = 0.0;
+
+ for(k=0;k<Cavity->nElem;k++)
+ {
+ newIndexes[k] = Cavity->elem[k];
+ AreaCavity += theMesh->Triangles[Cavity->elem[k]]->area;
+ }
+
+ for(k = Cavity->nElem; k < Edge1->nElem; k++)
+ {
+ newIndexes[k] = theMesh->nTriangles;
+ theMesh->nTriangles++;
+ }
+
+
+ // Etape 6 : Creation d'une structure contenant les pointeurs vers les nouveaux triangles
+ // allant remplir la cavite et dont un des sommets est le point (x,y). Les triangles d'aire
+ // negligeable ne seront pas ajoutes dans la triangulation (cela permet d'eviter d'eventuels
+ // problemes lorsque trois points sont presque alignes).
+
+ structTriangle **NewTriangles = malloc(sizeof(structTriangle*)*n);
+
+ int t1,t2,i,a=0,smallTri = 0;
+
+ for(i=0; i<n; i++)
+ {
+ t1 = Edge1->elem[i]; t2 = Edge2->elem[i];
+ if(t1>=0 && t2>=0)
+ {
+ structTriangle *Tin = theMesh->Triangles[t1];
+ structTriangle *Tout = theMesh->Triangles[t2];
+
+ int n0 = -2, n1 = -2, n2 = -2;
+ int j;
+ for(j=0;j<3;j++)
+ {
+ if (Tin->Nodes[0] == Tout->Nodes[j])
+ n0 = Tin->Nodes[0];
+
+ if (Tin->Nodes[1] == Tout->Nodes[j])
+ n1 = Tin->Nodes[1];
+
+ if (Tin->Nodes[2] == Tout->Nodes[j])
+ n2 = Tin->Nodes[2];
+ }
+
+
+ int p1,p2;
+ if(n0==-2) {p1 = n1; p2 = n2;}
+ else if(n1 == -2) {p1 = n2; p2 = n0;}
+ else if(n2 == -2) {p1 = n0; p2 = n1;}
+
+ NewTriangles[a] = TriangleNew(indexNode,p1,p2,theGrid);
+
+ NewTriangles[a]->Voisins[1] = Tout->index;
+
+ int t; for(t =0; t<3; t++)
+ {
+ if(Tout->Nodes[t] == p2)
+ {
+ Tout->Voisins[t] = newIndexes[a];
+ NewTriangles[a]->Segment[1] = Tout->Segment[t];
+ }
+ }
+ a++;
+ }
+ else
+ {
+ int nt, pt;
+ if(t1<0 && t2 >=0) {nt = t1; pt = t2;}
+ else if(t1>=0 && t2<0) {nt = t2; pt = t1;}
+ structTriangle* Tin = theMesh->Triangles[pt];
+ int n1 = -1, n2 = -1;
+ int f;
+ for(f=0;f<3;f++)
+ {
+ if (Tin->Voisins[f] == nt)
+ {
+ n1 = Tin->Nodes[f];
+ n2 = Tin->Nodes[(f+1)%3];
+ break;
+ }
+ }
+ structTriangle *myTriangle = TriangleNew(indexNode,n1,n2,theGrid);
+ myTriangle->Segment[1] = Tin->Segment[f];
+
+ double myArea = computeTriangleArea(myTriangle, theGrid);
+
+ // ne garder que les triangles suffisammment grands
+ if(myArea>=1e-7*fabs(AreaCavity))
+ {
+ NewTriangles[a] = myTriangle;
+ NewTriangles[a]->Voisins[1] = nt;
+ a++;
+ }
+ else
+ {
+ smallTri++;
+ NewTriangles[n-1] = myTriangle;
+ NewTriangles[n-1]->Voisins[1] = nt;
+ }
+
+ }
+ }
+
+
+ // Etape 7 : Mise a jour des voisins des nouveaux triangles et calcul de la surface
+ // totale des nouveaux triangles. Les noueveaux triangles sont egalement inseres
+ // dans l'arbre de tri.
+
+ double newArea = 0.0;
+ for(i=0;i<n;i++)
+ {
+ int v0;
+ for(v0=0;v0<n;v0++)
+ {
+ if(NewTriangles[i]->Nodes[1] == NewTriangles[v0]->Nodes[2])
+ {
+ NewTriangles[i]->Voisins[0] = newIndexes[v0];
+ NewTriangles[v0]->Voisins[2] = newIndexes[i];
+ break;
+ }
+ }
+ NewTriangles[i]->index = newIndexes[i];
+ if(i==n-1)
+ {
+ if(smallTri==0)
+ {
+ theMesh->Triangles[newIndexes[i]] = NewTriangles[i];
+ update_tree(NewTriangles[i],theMesh);
+ }
+ else
+ {
+ theMesh->nTriangles--;
+ }
+ }
+ else
+ {
+ theMesh->Triangles[newIndexes[i]] = NewTriangles[i];
+ update_tree(NewTriangles[i],theMesh);
+ }
+ newArea += NewTriangles[i]->area;
+
+ }
+
+ if(smallTri!=0)
+ {
+ theMesh->Triangles[NewTriangles[n-1]->Voisins[0]]->Voisins[2] = NewTriangles[n-1]->Voisins[1];
+ theMesh->Triangles[NewTriangles[n-1]->Voisins[2]]->Voisins[0] = NewTriangles[n-1]->Voisins[1];
+ theMesh->Triangles[NewTriangles[n-1]->Voisins[0]]->Segment[2] = NewTriangles[n-1]->Segment[1];
+ theMesh->Triangles[NewTriangles[n-1]->Voisins[2]]->Segment[0] = NewTriangles[n-1]->Segment[1];
+ free(NewTriangles[n-1]);
+ }
+
+
+ // Etape 8 : Verfier que la surface de la cavite est la meme avant et apres
+ // l'ajout des nouveaux triangles.
+
+ double deltaArea = fabs(newArea - AreaCavity);
+ if(deltaArea>(2e-7)*fabs(AreaCavity))
+ {
+ printf("Error: Cavity area changes \n");
+ printf("Old area = %le \n", AreaCavity);
+ printf("New area = %le \n", newArea);
+ exit(0);
+ }
+
+
+ // Etape 9 : Liberation de la memoire.
+
+ free(newIndexes);
+ freeVector(Edge1);
+ freeVector(Edge2);
+ freeVector(Cavity);
+ free(NewTriangles);
+
+}
+
+/*
+OBJECTIFS:
+Repertorier les aretes du contour de la cavite. L'alorithme parcourt toutes
+les aretes de tous les triangles de la cavite. Si une arete se trouve deja
+dans la liste des aretes du contour de la cavite, alors cette arete est supprimee
+de la liste (cad la derniere arete de la liste est mise a la place de cette arete)
+sinon, l'arete est ajoutee dans la liste. Dans les deux cas, le nombre d'aretes
+de la liste finale est adapte.
+Une arete est definie par les deux triangles adjacents a cette arete.
+
+INPUTS:
+-> Edge1 : vecteur vide de taille 3 x nombre de triangles de la cavite.
+-> Edge2 : vecteur vide de taille 3 x nombre de triangles de la cavite.
+-> Cavity : structure contenant des pointeurs vers tous les triangles de la cavite.
+-> theTriangulation : structure contenant des pointeurs vers tous les triangles
+ de la triangulation courante.
+
+OUTPUTS:
+-> Mise à jour de Edge1 : contient des pointeurs vers les triangles à gauche des
+ aretes du contour de la cavite.
+-> Mise à jour de Edge2 : contient des pointeurs vers les triangles à droite des
+ aretes du contour de la cavite.
+*/
+
+void GetEdges(structVector *Edge1, structVector *Edge2, structVector *Cavity, structTriangulation *theTriangulation)
+{
+ int i,v,j,n,m,toAdd;
+
+ // Etape 1 : Parcourir tous les triangles et toutes les aretes de la cavite.
+
+ for(j=0;j<Cavity->nElem;j++)
+ {
+ i = Cavity->elem[j];
+
+ for(n=0;n<3;n++)
+ {
+ toAdd = 1;
+ v = theTriangulation->Triangles[i]->Voisins[n];
+
+ // Etape 2 : Si l'arete courante se trouve deja dans la liste alors
+ // supprimer cette arete de la liste (ecrire la derniere arete de la
+ // liste a sa place).
+
+ for(m=0; m<Edge1->nElem; m++)
+ {
+ if(i == Edge2->elem[m] && v == Edge1->elem[m])
+ {
+ Edge1->elem[m] = Edge1->elem[Edge1->nElem-1];
+ Edge2->elem[m] = Edge2->elem[Edge2->nElem-1];
+ Edge1->nElem --; Edge2->nElem--;
+ toAdd = 0; break;
+ }
+ }
+
+ // Etape 3 : Si l'arete courante ne se trouve pas encore dans la liste
+ // alors elle peut y etre ajoutee.
+
+ if(toAdd)
+ {
+ Edge1->elem[Edge1->nElem] = i;
+ Edge2->elem[Edge2->nElem] = v;
+ Edge1->nElem ++; Edge2->nElem++;
+ }
+ }
+ }
+}
+
+
+/*
+OBJECTIFS:
+Determiner les triangles faisant partie de la cavite associe au point (x,y).
+L'algorithme part du triangle dans lequel est contenu le point (x,y). Il parcourt
+ensuite les voisins de ce triangle de proche en proche et ajoute dans la cavité
+tous les triangles dont le cercle circonscrit contient le poinr (x,y).
+L'algorithme veille a ne pas tester deux fois le meme triangle grace a la variable
+"flagged" qui vaut 1 si un triangle a deja ete teste.
+
+INPUTS:
+-> Cavity : vecteur initialement vide (la taille est inconnue a priori) et sera
+ fixee arbitrairement dans la fonction addNode.
+-> theTriangulation : structure contenant des pointeurs vers tous les triangles
+ de la triangulation courante.
+-> T : pointeur vers le triangle a tester.
+-> x et y : coordonnees du point associe a la cavite.
+-> theGrid: contient les coordonnees de l'ensemble des points fournis.
+
+OUTPUTS:
+-> La structure Cavity est mise a jour et contient des pointeurs vers tous les
+ triangles de la cavite.
+*/
+
+void buildCavity(structVector *Cavity, structTriangulation *theTriangulation, structTriangle *T, double x, double y, structGrid *theGrid)
+{
+ if(inCircle(x,y,T,theGrid))
+ {
+ // Etape 1 : Si le cercle circonscrit du triangle courant contient le point(x,y)
+ // alors ce triangle peut etre ajoute dans la cavite.
+
+ Cavity->elem[Cavity->nElem] = T->index;
+ Cavity->nElem ++;
+
+ // Etape 2 : Indiquer a l'algorithme de ne plus tester ce triangle.
+
+ T->flagged = 1;
+
+ // Etape 3 : Tester les triangles adjacents a chaque triangle ajoute dans la cavite.
+
+ int i,indexVoisin;
+ structTriangle *V;
+ for(i=0;i<3;i++)
+ {
+ indexVoisin = T->Voisins[i];
+ if (indexVoisin > -1)
+ {
+ V = theTriangulation->Triangles[indexVoisin];
+ if (!V->flagged)
+ buildCavity(Cavity, theTriangulation, V,x,y,theGrid);
+ }
+ }
+ }
+}
+
+
+/*
+OBJECTIFS:
+Lorsqu'un noeud a ete ajoute dans la triangulation, il est necesssaire de
+verifier que ce noeud n'empiete sur aucun segment. Au lieu de tester l'ensemble
+des segments, nous allons appliquer une heuristique.
+Cette methode se base sur le principe que tous les segments empietes par un noeud
+seront compris dans la cavite associee a ce point en question.
+La cavite associe au noeud est calculee. Toutes les aretes des triangles de la cavite
+seront testees. Parmis les aretes, celles qui correspondent a des segments dont le
+cercle diametral contient le noeud considere seront memorisees.
+
+
+INPUTS:
+-> encSegs : vecteur vide
+-> theMesh : structure contenant des pointeurs vers tous les triangles
+ de la triangulation courante.
+-> W : pointeur vers un triangle susceptible d'etre proche de la cavite.
+-> x et y : coordonnees du point associe a la cavite.
+-> theGrid: contient les coordonnees de l'ensemble des points fournis.
+-> seg1 et seg2 : liste des extremites des segments.
+
+OUTPUTS:
+-> encSegs : contient les indices des segments qui sont empietes par le noeud (x,y)
+*/
+void CavityTest(structVector *encSegs, structTriangulation *theMesh, structTriangle *W, double x, double y, structGrid *theGrid,structVector *seg1, structVector *seg2)
+{
+ // Etape 1 : Reperer tous les triangles de la cavite associee au point (x,y)
+
+ structTriangle *Tri = FindFirstTri(W, x, y, theGrid, theMesh,-1,1);
+ structVector *Cavity = newVector(1000);
+ buildCavity(Cavity,theMesh,Tri,x,y,theGrid);
+ int i = 0, j = 0;
+ for(i = 0; i<Cavity->nElem; i++)
+ {
+ structTriangle *T = theMesh->Triangles[Cavity->elem[i]];
+ T->flagged = 0;
+ }
+
+
+ // Etape 2 : Tester toutes les aretes de tous les triangles de la cavite.
+ // Reperer les aretes qui sont egalement des segments. Tester si le noeud
+ // est contenu dans le cercle diametral associe a ce segment.
+
+ for(i = 0; i<Cavity->nElem; i++)
+ {
+ structTriangle *T = theMesh->Triangles[Cavity->elem[i]];
+ for(j = 0; j<3; j++)
+ {
+ int ind = T->Segment[j];
+ if(ind != -1)
+ {
+ if(checkSeg(seg1,seg2, theGrid, x, y, ind))
+ {
+ int alreadyDone = 0;
+ int n;
+ for(n=0;n<encSegs->nElem;n++)
+ {
+ if(ind == encSegs->elem[n])
+ alreadyDone = 1;
+ }
+ if(!alreadyDone)
+ {
+ encSegs->elem[encSegs->nElem] = ind;
+ encSegs->nElem++;
+ }
+ }
+ }
+ }
+ }
+ freeVector(Cavity);
+}
+
+
+/*
+OBJECTIFS:
+Le but de cette fonction est de localiser le point (x,y). Cette localisation peut
+se faire suivant deux critères:
+-> Soit trouver exactement dans quel triangle se trouve le point (x,y) : circle = 0
+-> Soit trouver un triangle dont le cercle circonscrit contient le point (x,y): circle = 1
+La recherche se base sur la methode de la droite. Un vecteur est trace entre un
+triangle de depart et le point (x,y). Pour savoir quel voisin du triangle doit etre
+teste, on calcule deux produits vectoriels:
+-> produit vectoriel entre un vecteur reliant le barycentre du triangle a un de ses sommets
+ (A) et le vecteur de recherche.
+-> produit vectoriel entre le vecteur de recherche et le vecteur reliant le barycentre du
+ triangle a un autre de ses sommets (B).
+Si les deux produits vectoriels sont positifs, alors l'algorithme peut tester le triangle
+adjacent a l'arete A-B.
+
+
+INPUTS:
+-> T : pointeur vers un triangle susceptible d'etre proche du point a localiser.
+-> x et y : coordonnees du point a localiser.
+-> theGrid: contient les coordonnees de l'ensemble des points fournis.
+-> theMesh : structure contenant des pointeurs vers tous les triangles
+ de la triangulation courante.
+-> prevTri : indice du precedent triangle visite.
+-> circle : critere de recherche: cerle circonscrit (1) ou triangle (0).
+
+
+OUTPUTS:
+-> Renvoie un pointeur vers un triangle "contenant" le poinr (x,y)
+*/
+structTriangle* FindFirstTri(structTriangle *T, double x, double y, structGrid *theGrid, structTriangulation *theMesh, int prevTri,int circle)
+{
+ // Etape 1 : verifier si le triangle a ete trouve
+
+ if(circle == 1)
+ {
+ if(inCircle(x,y,T,theGrid))
+ {
+ return T;
+ }
+ }
+ else if (circle==0)
+ {
+ if(isInTriangle(theGrid,T,x,y))
+ {
+ return T;
+ }
+ }
+
+
+ // Etape 2 : si le triangle n'a pas ete trouve, chercher
+ // le voisin de ce triangle qui devra etre teste.
+
+
+ // Etape 2a : calcul du barycentre du triangle courant.
+ double nodesX[3];
+ double nodesY[3];
+
+ nodesX[0] = theGrid->X[T->Nodes[0]-1];
+ nodesY[0] = theGrid->Y[T->Nodes[0]-1];
+ nodesX[1] = theGrid->X[T->Nodes[1]-1];
+ nodesY[1] = theGrid->Y[T->Nodes[1]-1];
+ nodesX[2] = theGrid->X[T->Nodes[2]-1];
+ nodesY[2] = theGrid->Y[T->Nodes[2]-1];
+
+ double nx = (nodesX[0]+nodesX[1]+nodesX[2])/3.0;
+ double ny = (nodesY[0]+nodesY[1]+nodesY[2])/3.0;
+
+ // Etape 2b: calcul des produits vectoriels.
+ double ax = x - nx;
+ double ay = y - ny;
+
+ int i;
+ for(i = 0; i<3; i++)
+ {
+ double v1x = nodesX[i] - nx;
+ double v1y = nodesY[i] - ny;
+ double v2x = nodesX[(i+1)%3] - nx;
+ double v2y = nodesY[(i+1)%3] - ny;
+
+ if(ax*v1y-ay*v1x > -1e-15*(ax*ax+ay*ay) && ax*v2y-ay*v2x < 1e-15*(ax*ax+ay*ay))
+ {
+ if (T->Voisins[i] < 0 || T->Voisins[i] == prevTri)
+ {
+ printf("Error findFirstTri \n");
+ printf("voisin = %d previous = %d \n", T->Voisins[i], prevTri);
+ exit(0);
+ }
+ else
+ {
+ return FindFirstTri(theMesh->Triangles[T->Voisins[i]],x,y,theGrid,theMesh,T->index,circle);
+ }
+ }
+ }
+ printf("No direction found \n");
+ exit(0);
+}
+
+
diff --git a/addTriangle.c b/addTriangle.c
new file mode 100644
index 0000000..0ee4ca9
--- /dev/null
+++ b/addTriangle.c
@@ -0,0 +1,75 @@
+#include "header.h"
+
+/*
+OBJECTIFS:
+addTriangle est une fonction permettant d'inserer un triangle dans un arbre qui les
+trie selon leur qualite. La structure d'arbre permet également de creer une liste
+chainee demarrant au plus mauvais triangle, reference par begin.
+
+INPUTS:
+-> theTriangle, le triangle a inserer
+-> node, l'arbre dans lequel on insere theTriangle
+-> theMesh, la triangulation actuelle
+-> les noeuds les plus a droite et a gauche du noeud a inserer a une iteration donnee
+
+OUTPUTS:
+-> l'arbre est mis a jour
+-> le pointeur begin de theMesh est mis a jour
+*/
+
+void addTriangle(structTriangle *theTriangle,
+ structQualityTree *node,\
+ structTriangulation *theMesh,\
+ structQualityTree *previousL,\
+ structQualityTree *previousR)
+{
+ // Si le noeud ne reference pas encore un triangle, on ajoute le triangle dans l'arbre
+ if (node->quality < 0.0)
+ {
+ node->quality = theTriangle->quality;
+ node->triangle = theTriangle;
+ node->left = QualityTreeNew();
+ node->right = QualityTreeNew();
+
+ // On insere l'element dans la liste chainee
+ if(previousL != NULL)
+ previousL->next = node;
+
+ node->next = previousR;
+
+ // On met begin a jour
+ if (theMesh->begin->quality <= node->quality)
+ {
+ theMesh->begin = node;
+ }
+ }
+
+ // Si le noeud reference deja un triangle, on descend dans l'arbre
+ else
+ {
+ // Si le triangle a une qualite plus mauvaise, on descend a gauche
+ if (theTriangle->quality >= node->quality)
+ {
+ previousR = node;
+ addTriangle(theTriangle,node->left,theMesh,previousL,previousR);
+ }
+
+ // Si le triangle a une qualite meilleure, on descend a gauche
+ else
+ {
+ previousL = node;
+ addTriangle(theTriangle,node->right,theMesh,previousL,previousR);
+ }
+ }
+}
+
+
+/*
+OBJECTIFS:
+Appel initial (avant recursion) de la fonction addTriangle (ci-dessus)
+*/
+
+void update_tree(structTriangle *theTriangle, structTriangulation * theTriangulation)
+{
+ addTriangle(theTriangle,theTriangulation->tree,theTriangulation,NULL,NULL);
+}
diff --git a/checkAndSplitEncroachedSegments.c b/checkAndSplitEncroachedSegments.c
new file mode 100644
index 0000000..46bc43d
--- /dev/null
+++ b/checkAndSplitEncroachedSegments.c
@@ -0,0 +1,281 @@
+#include "header.h"
+/*
+OBJECTIFS:
+Cette fonction permet de s'assurer que tous les segments sont dans triangulation.
+Pour cela, pour chaque segement:
+1) on verifie si un noeud est dans le cercle diametral de ce segment. Au lieu de
+ parcourir toute la liste de noeuds, on fait appel efficacement a la fonction "toSplit".
+2) Si le segment est effictement empiete par un noeud, alors on calcule les
+ coordonnnes du point de division du segment a l'aide de la fonction "setSplittingPoint".
+3) Le problème est que ce nouveau noeud pourrait empieter sur d'autres segments. Pour verifier
+ si ce noeud empiete sur un segment, au lieu de parcourir tous les segments, on utilise la
+ fonction "cavityTest".
+4) Enfin, le segment est divise en deux sous segements. Les deux sous segments sont ajoutes
+ a la liste des segments et leur validite sera aussi testee.
+
+*/
+void checkOneSeg(int i, int indexCloseTri, structVector *seg1,structVector *seg2,structGrid *theGrid,structTriangulation *theMesh)
+{
+ int n[2]; n[0] = seg1->elem[i]; n[1] = seg2->elem[i];
+
+ double x1 = theGrid->X[n[0]-1]-1e-8*Sign(theGrid->X[n[0]-1]);
+ double y1 = theGrid->Y[n[0]-1]-1e-8*Sign(theGrid->Y[n[0]-1]);
+
+ structTriangle *T = FindFirstTri(theMesh->Triangles[indexCloseTri], x1, y1, theGrid, theMesh,-1,1);
+ int indCloseTri = T->index;
+
+ structVector *EncSegs = newVector(1000);
+
+ int j;
+ for(j = 0; j<2; j++)
+ {
+ int node = n[j];
+ double x = theGrid->X[node-1];
+ double y = theGrid->Y[node-1];
+
+ double shift1 = fabs(rand()/(double)RAND_MAX)*1e-8*Sign(x);
+ double shift2 = fabs(rand()/(double)RAND_MAX)*1e-8*Sign(y);
+ T = FindFirstTri(theMesh->Triangles[indCloseTri], x-shift1, y-shift2, theGrid, theMesh,-1,0);
+
+ int z; int inTri = 0;
+ for(z = 0; z<3; z++)
+ {
+ if(T->Nodes[z] == node)
+ inTri = 1;
+ }
+
+ if(inTri == 0)
+ {
+ CavityTest(EncSegs, theMesh, theMesh->Triangles[indCloseTri], x, y, theGrid, seg1, seg2);
+ addNode(theMesh->Triangles[indCloseTri], node,theMesh, theGrid,seg1,seg2);
+ }
+ }
+
+ int N = seg1->nElem; int k;
+ for(k=i; k<seg1->nElem;)
+ {
+ if(toSplit(theMesh->Triangles[indCloseTri], seg1, seg2, k, theMesh, theGrid) != -1)
+ {
+ double x,y;
+ setSplittingPoint(theGrid, seg1, seg2, k, &x, &y,1);
+ CavityTest(EncSegs,theMesh, theMesh->Triangles[indCloseTri], x, y, theGrid, seg1, seg2);
+ splitSeg(theMesh->Triangles[indCloseTri],theMesh,theGrid,seg1,seg2,k,1);
+ }
+ else
+ {
+ updateSegments(theMesh->Triangles[indCloseTri],seg1, seg2, k, theMesh,theGrid);
+ if (k == i) k = N; else k++;
+ }
+ }
+
+ for(k=0; k<EncSegs->nElem; k++)
+ {
+ double x,y;
+ setSplittingPoint(theGrid, seg1, seg2, EncSegs->elem[k], &x, &y,1);
+ CavityTest(EncSegs,theMesh, theMesh->Triangles[i], x, y, theGrid, seg1, seg2);
+ splitSeg(theMesh->Triangles[indCloseTri],theMesh,theGrid,seg1,seg2,EncSegs->elem[k],1);
+
+ int y1 = toSplit(theMesh->Triangles[i], seg1, seg2, EncSegs->elem[k], theMesh, theGrid);
+ if(y1 == -1)
+ {
+ updateSegments(theMesh->Triangles[i],seg1, seg2, EncSegs->elem[k], theMesh,theGrid);
+ }
+ else
+ {
+ EncSegs->elem[EncSegs->nElem] = EncSegs->elem[k]; EncSegs->nElem ++;
+ }
+
+ int y2 = toSplit(theMesh->Triangles[i], seg1, seg2, seg2->nElem-1, theMesh, theGrid);
+ if(y2 == -1)
+ {
+ updateSegments(theMesh->Triangles[i],seg1, seg2, seg2->nElem-1, theMesh,theGrid);
+ }
+ else
+ {
+ EncSegs->elem[EncSegs->nElem] = seg2->nElem-1; EncSegs->nElem ++;
+ }
+ EncSegs->elem[k] = -1;
+ }
+
+ freeVector(EncSegs);
+}
+
+/*
+OBJECTIFS:
+L'objectif de cette fonction est de parcourir l'ensemble des segments, de verifier
+si un de ces segments est empietes par un noeud et le cas echeant, de diviser ce
+segment.
+*/
+int checkAndSplitEncroachedSegments(structVector *seg1,structVector *seg2,structGrid *theGrid,structTriangulation *theMesh, structVector *segMod)
+{
+ int N = seg1->nElem;
+ int i;
+
+
+ for (i=0; i<N; ++i)
+ {
+ int nTri = theMesh->nTriangles;
+ checkOneSeg(i, nTri-1, seg1, seg2, theGrid, theMesh);
+ }
+}
+
+/*
+OBJECTIFS:
+Pour un segment donne, on veut verifier rapidement si ce segment est
+empiete par un noeud, sans parcourir toute la liste de noeuds. Pour cela,
+on teste les deux triangles adjacents a ce segment.
+Pour chacun de ces deux triangles:
+-> si le sommet opppose au segment empiete sur ce segment, alors on sait
+ qu'il faudra couper le segment et l'objectif de la fonction est atteint.
+-> Par contre, si le sommet oppose au segment n'empiete pas sur ce segment,
+ alors on peut cerifier qu'aucun autre noeud n'empiete sur ce segment.
+ En effet, dans ce cas, le cercle diametral sera inclus dans le cercle
+ circonscrit au triangle. Or, puisque la triangulation de Delaunay est
+ maintenue tout au long l'algorithme, il n'y aura aucun point dans le
+ cercle circonscrit, et donc aucun point dans le cercle diametral.
+*/
+int toSplit(structTriangle *W, structVector *seg1, structVector *seg2, int n, structTriangulation* theMesh, structGrid* theGrid)
+{
+ int n1 = seg1->elem[n];
+ int n2 = seg2->elem[n];
+
+ double n1x = theGrid->X[n1-1];
+ double n1y = theGrid->Y[n1-1];
+ double n2x = theGrid->X[n2-1];
+ double n2y = theGrid->Y[n2-1];
+
+ double normaleX = -(n2y-n1y);
+ double normaleY = n2x-n1x;
+
+ double mx = 0.5*(n1x+n2x);
+ double my = 0.5*(n1y+n2y);
+
+ if(normaleX*mx+normaleY*my > 0)
+ {
+ normaleX = -normaleX;
+ normaleY = -normaleY;
+ }
+ mx = mx+normaleX*1e-8;
+ my = my+normaleY*1e-8;
+
+ structTriangle *T = FindFirstTri(W, mx, my, theGrid,theMesh, -1, 0);
+
+
+ int i;
+ for(i=0; i<3; i++)
+ {
+ if(n1 == T->Nodes[i] && n2 == T->Nodes[(i+1)%3])
+ {
+ break;
+ }
+ if(n2 == T->Nodes[i] && n1 == T->Nodes[(i+1)%3])
+ {
+ break;
+ }
+ if(i==2)
+ {
+ return T->index;
+ }
+ }
+
+ int j = T->Nodes[(i+2)%3];
+
+ double x = theGrid->X[j-1];
+ double y = theGrid->Y[j-1];
+
+ if(inDiametralCircle(n1x, n1y, n2x, n2y, x, y))
+ return T->index;
+
+ int k = T->Voisins[i];
+
+ if (k<0)
+ {
+ return -1;
+ }
+
+ structTriangle *V = theMesh->Triangles[k];
+
+ for(i=0; i<3; i++)
+ {
+ if(V->Voisins[i] == T->index)
+ {
+ j = V->Nodes[(i+2)%3];
+
+ x = theGrid->X[j-1];
+ y = theGrid->Y[j-1];
+
+ if(inDiametralCircle(n1x, n1y, n2x, n2y, x, y))
+ return T->index;
+ else
+ {
+ return -1;
+ }
+ }
+ }
+}
+
+/*
+OBJECTIFS:
+Cette fonction permet de mettre a jour les segments. Lorsqu'un segment a ete divise,
+ses deux sous segments ont ete ajoutes dans la liste des segments. Cependant, il est
+necessaire de parcourir les triangles adjacents a ces segments, afin d'indiquer quelles
+aretes de ces triangles sont maintenant des segments. En bref, il s'agit de mettre
+a jour le vecteur "Segment" des structures "structTriangles".
+*/
+void updateSegments(structTriangle *W, structVector *seg1, structVector *seg2, int n, structTriangulation* theMesh, structGrid* theGrid)
+{
+ int n1 = seg1->elem[n];
+ int n2 = seg2->elem[n];
+
+ double n1x = theGrid->X[n1-1];
+ double n1y = theGrid->Y[n1-1];
+ double n2x = theGrid->X[n2-1];
+ double n2y = theGrid->Y[n2-1];
+
+ double normaleX = -(n2y-n1y);
+ double normaleY = n2x-n1x;
+
+ double mx = 0.5*(n1x+n2x);
+ double my = 0.5*(n1y+n2y);
+
+ if(normaleX*mx+normaleY*my > 0)
+ {
+ normaleX = -normaleX;
+ normaleY = -normaleY;
+ }
+
+ mx = mx+normaleX*1e-8;
+ my = my+normaleY*1e-8;
+
+
+ structTriangle *T = FindFirstTri(W, mx, my, theGrid,theMesh, -1, 0);
+
+ int i;
+ for(i=0; i<3; i++)
+ {
+ if(n1 == T->Nodes[i] && n2 == T->Nodes[(i+1)%3])
+ {
+ break;
+ }
+ if(n2 == T->Nodes[i] && n1 == T->Nodes[(i+1)%3])
+ {
+ break;
+ }
+ }
+
+ T->Segment[i] = n;
+
+ if(T->Voisins[i] >= 0)
+ {
+ structTriangle *V = theMesh->Triangles[T->Voisins[i]];
+
+ for(i=0; i<3; i++)
+ {
+ if(V->Voisins[i] == T->index)
+ {
+ V->Segment[i] = n;
+ break;
+ }
+ }
+ }
+}
diff --git a/checkAndSplitSkinnyTriangle.c b/checkAndSplitSkinnyTriangle.c
new file mode 100644
index 0000000..93cd96d
--- /dev/null
+++ b/checkAndSplitSkinnyTriangle.c
@@ -0,0 +1,164 @@
+#include "header.h"
+
+/*
+OBJECTIFS:
+Cette fonction permet de supprimer un triangle de mauvaise qualite. Si le centre du cercle circonscrit
+n'empiete sur aucun segment, alors le centre du cercle circonscrit est ajoute dans la triangulation.
+Sinon, pour chaque segment empiete:
+1) on calcule les coordonnees du point qui divisera le segment
+2) on verifie que ce nouveau point n'empiete sur aucun autre segment. Si jamais un nouveau segment est
+ empiete, alors ce segment est ajoute a la liste des segments empietes.
+3) on divise le segment
+4) on met a jour la liste des segments.
+*/
+void checkAndSplitSkinnyTriangle(structVector *seg1,structVector *seg2,structGrid *theGrid,structTriangulation *theMesh, int i)
+{
+ double px,py;
+ computeCircumcenter(theMesh->Triangles[i],theGrid,&px,&py);
+
+ structVector *EncSegs = newVector(300);
+ CavityTest(EncSegs, theMesh, theMesh->Triangles[i], px, py, theGrid, seg1, seg2);
+
+ if (EncSegs->nElem == 0)
+ {
+ splitTri(theMesh->Triangles[i],theMesh, theGrid, px,py,seg1,seg2);
+ }
+
+ int k; double x,y;
+ int N = EncSegs->nElem;
+
+ for(k=0; k<EncSegs->nElem; k++)
+ {
+ setSplittingPoint(theGrid, seg1, seg2, EncSegs->elem[k], &x, &y,k>=N);
+ CavityTest(EncSegs,theMesh, theMesh->Triangles[i], x, y, theGrid, seg1, seg2);
+
+ splitSeg(theMesh->Triangles[i],theMesh,theGrid,seg1,seg2,EncSegs->elem[k],k>=N);
+
+ int y1 = toSplit(theMesh->Triangles[i], seg1, seg2, EncSegs->elem[k], theMesh, theGrid);
+ if(y1 == -1)
+ {
+ updateSegments(theMesh->Triangles[i],seg1, seg2, EncSegs->elem[k], theMesh,theGrid);
+ }
+ else
+ {
+ EncSegs->elem[EncSegs->nElem] = EncSegs->elem[k]; EncSegs->nElem ++;
+ }
+
+ int y2 = toSplit(theMesh->Triangles[i], seg1, seg2, seg2->nElem-1, theMesh, theGrid);
+ if(y2 == -1)
+ {
+ updateSegments(theMesh->Triangles[i],seg1, seg2, seg2->nElem-1, theMesh,theGrid);
+ }
+ else
+ {
+ EncSegs->elem[EncSegs->nElem] = seg2->nElem-1; EncSegs->nElem ++;
+ }
+ EncSegs->elem[k] = -1;
+ }
+ free(EncSegs);
+}
+
+/*
+OBJECTIFS:
+Cette fonction permet de traiter le cas des petits angles presents dans le fichier d'input.
+Lorsqu'un petit angle est compris entre deux segments et que ces segments sont de longueur
+minimum, alors le petit triangle est supprime. Pour cela, l'arete opposee au petit angle
+devient un segment et les deux autres segments sont supprimes. On verifie ensuite que le nouveau
+segment n'est empiete par aucun point.
+*/
+int canBeSplit(structTriangle *T, structVector *seg1, structVector *seg2, structGrid *theGrid, double minAngle, structTriangulation *theMesh)
+{
+ double cosMax = cos(minAngle);
+
+ int j;
+ for(j=0; j<3; j++)
+ {
+ if(T->Segment[j] >= 0 && T->Segment[(j+1)%3] >= 0)
+ {
+ int s1 = T->Segment[j];
+ int s2 = T->Segment[(j+1)%3];
+
+ int m,n1,n2;
+ if(seg1->elem[s1] == seg1->elem[s2])
+ {
+ m = seg1->elem[s1];
+ n1 = seg2->elem[s1];
+ n2 = seg2->elem[s2];
+ }
+ else if(seg2->elem[s1] == seg2->elem[s2])
+ {
+ m = seg2->elem[s1];
+ n1 = seg1->elem[s1];
+ n2 = seg1->elem[s2];
+ }
+ else if(seg1->elem[s1] == seg2->elem[s2])
+ {
+ m = seg1->elem[s1];
+ n1 = seg2->elem[s1];
+ n2 = seg1->elem[s2];
+ }
+ else if(seg2->elem[s1] == seg1->elem[s2])
+ {
+ m = seg2->elem[s1];
+ n1 = seg1->elem[s1];
+ n2 = seg2->elem[s2];
+ }
+ else
+ {
+ printf("Error in corner lopping \n");
+ exit(0);
+ }
+ double n1x = theGrid->X[n1-1];
+ double n1y = theGrid->Y[n1-1];
+ double n2x = theGrid->X[n2-1];
+ double n2y = theGrid->Y[n2-1];
+ double mx = theGrid->X[m-1];
+ double my = theGrid->Y[m-1];
+
+ double s1x = n1x - mx;
+ double s1y = n1y - my;
+ double s2x = n2x - mx;
+ double s2y = n2y - my;
+
+ double length1 = sqrt(s1x*s1x+s1y*s1y);
+ double length2 = sqrt(s2x*s2x+s2y*s2y);
+
+ double cosine = (s1x*s2x+s1y*s2y)/(length1*length2);
+
+ if(cosine >= cosMax && (length1<=0.5*seg1->lengthMinSeg || length2<=0.5*seg1->lengthMinSeg))
+ {
+ for(j = 0; j<3; j++)
+ {
+ T->Segment[j] = -1;
+ int v = T->Voisins[j];
+
+ if(v>=0)
+ {
+ structTriangle *V = theMesh->Triangles[v];
+
+ int z;
+
+ for (z = 0; z<3; z++)
+ {
+ if(V->Voisins[z] == T->index)
+ {
+ V->Segment[z] = -1;
+ break;
+ }
+ }
+ }
+ }
+
+ seg1->elem[s1] = n1;
+ seg2->elem[s1] = n2;
+
+
+ checkOneSeg(s1, T->index, seg1, seg2, theGrid, theMesh);
+ //checkOneSeg(s2, T->index, seg1, seg2, theGrid, theMesh);
+
+ return 1;
+ }
+ }
+ }
+ return 1;
+}
diff --git a/freeStructures.c b/freeStructures.c
new file mode 100644
index 0000000..5cdd48b
--- /dev/null
+++ b/freeStructures.c
@@ -0,0 +1,67 @@
+#include "header.h"
+
+/*
+OBJECTIF:
+Liberer la memoire allouee pour une structure de type vecteur.
+*/
+
+void freeVector(structVector *theVector)
+{
+ free(theVector->elem);
+ theVector->elem = NULL;
+ free(theVector);
+ theVector = NULL;
+}
+
+
+/*
+OBJECTIF:
+Liberer la memoire allouee pour une stcuture de type triangulation.
+Seuls les pointeurs vers les triangles sont supprimes. Les triangles
+en eux-memes seront liberes lorsque l'arbre sera libere.
+*/
+
+void freeTriangulation(structTriangulation *theMesh)
+{
+ free(theMesh->Triangles);
+ theMesh->Triangles = NULL;
+ free(theMesh);
+ theMesh = NULL;
+}
+
+
+/*
+OBJECTIF:
+Liberer la memoire allouee pour une structure de type Grid.
+*/
+
+void freeGrid(structGrid *theGrid)
+{
+ free(theGrid->X);
+ free(theGrid->Y);
+ free(theGrid->Z);
+ theGrid->X = NULL;
+ theGrid->Y = NULL;
+ theGrid->Z = NULL;
+ free(theGrid);
+}
+
+/*
+OBJECTIF:
+Liberer l'arbre qui trie les triangles suivant leur qualite.
+La liberation de la memoire se fait de maniere recursive a partir
+des feuilles de l'arbre.
+*/
+void freeQualityTree(structQualityTree *tree)
+{
+ if(tree->quality > 0.0)
+ {
+ freeQualityTree(tree->left);
+ freeQualityTree(tree->right);
+ if(tree->triangle!=NULL)
+ free(tree->triangle);
+
+ }
+ free(tree);
+ tree = NULL;
+}
diff --git a/header.h b/header.h
new file mode 100644
index 0000000..794e0c4
--- /dev/null
+++ b/header.h
@@ -0,0 +1,130 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#define _GLFEM_H_
+#include <GL/glfw.h>
+#include <time.h>
+
+
+# define Min(a,b) ((a < b) ? a : b )
+# define Max(a,b) ((a > b) ? a : b )
+# define Abs(a) ((a > 0) ? a : -(a))
+# define Sign(a) ((a > 0) ? 1 : -1)
+
+
+typedef struct structTriangle structTriangle;
+struct structTriangle
+{
+ int Nodes[3];
+ int Voisins[3];
+ int index;
+ int flagged;
+ double area;
+ double quality;
+ int Segment[3];
+ int toBeDeleted;
+};
+
+typedef struct structQualityTree structQualityTree;
+struct structQualityTree{
+ structQualityTree *left;
+ structQualityTree *right;
+ structQualityTree *next;
+ double quality;
+ structTriangle *triangle;
+};
+
+typedef struct {
+ double *X;
+ double *Y;
+ double *Z;
+ int nNode;
+ int nInputVertex;
+ int size;
+} structGrid;
+
+typedef struct {
+ structTriangle **Triangles;
+ int nTriangles;
+ structQualityTree * tree;
+ structQualityTree * begin;
+ int size;
+} structTriangulation;
+
+typedef struct {
+ int *elem;
+ int nElem;
+ int size;
+ double lengthMinSeg;
+} structVector;
+
+
+
+// Initialisation des structures
+structVector *newVector(int size);
+structTriangulation *TriangulationNew(structGrid *theGrid);
+structTriangle *TriangleNew(int n1,int n2, int n3, structGrid *theGrid);
+structQualityTree * QualityTreeNew();
+
+// Lecture du fichier d'input
+void readInputFile(structGrid *theGrid, const char *filename);
+void shuffle(structGrid *theGrid);
+
+// Traitement des 4 points externes
+void addExteriorNodes(structGrid *theGrid,structVector *seg1,structVector *seg2,structTriangulation *theTriangulation);
+
+// Ajout de points (methode de la cavite)
+void addNode(structTriangle *theTriangle, int indexNode, structTriangulation *theMesh, structGrid *theGrid, structVector *seg1, structVector *seg2);
+void GetEdges(structVector *EdgeN1, structVector *EdgeN2, structVector *Cavity, structTriangulation *theMesh);
+void buildCavity(structVector *Cavity, structTriangulation *theMesh, structTriangle *T,double x, double y,structGrid *theGrid);
+
+// Tester si un point est dans un triangle ou dans un cercle
+int checkSeg(structVector *seg1,structVector *seg2,structGrid *theGrid, double x, double y, int ind);
+int isInTriangle(structGrid *theGrid, structTriangle *theTriangle, double x, double y);
+int sameSide(double x1, double y1, double x2, double y2, double xp, double yp, double xq, double yq);
+int inCircle(double dx, double dy, structTriangle *T,structGrid* theGrid);
+double computeTriangleArea(structTriangle *theTriangle,structGrid *theGrid);
+int inDiametralCircle(double ax, double ay, double bx, double by, double cx, double cy);
+
+// Interface graphique
+void DrawCircle(structTriangle *Triangle, double *X, double *Y);
+void showIt(structGrid *theGrid, structTriangulation *theTriangulation,structVector *seg1,structVector *seg2,int Dimension2, int withCircles, double angle);
+void Shut_Down(int return_code);
+void Draw_Triangulation(int nNodes, double *x, double *y, double *z, structTriangulation *T, int nTriang,int Dimension2, int withCircles,structVector *seg1,structVector *seg2, double angle);
+
+double min(double *x, int n);
+double max(double *x, int n);
+float ReshapeAndShift(structGrid *theGrid);
+
+// Liberation de la memoire allouee
+void freeQualityTree(structQualityTree *tree);
+void freeVector(structVector *theVector);
+void freeGrid(structGrid *theGrid);
+void freeTriangulation(structTriangulation *theMesh);
+
+// Fonction permettant d'afficher le README à la console
+void printManPage(const char *);
+
+// Nouvelles fonctions
+int canBeSplit(structTriangle *T, structVector *seg1, structVector *seg2, structGrid *theGrid, double minAngle, structTriangulation *theMesh);
+void setSplittingPoint(structGrid *theGrid, structVector *seg1, structVector *seg2, int i, double *x, double *y, int a);
+void CavityTest(structVector *encSegs, structTriangulation *theMesh, structTriangle *W, double x, double y, structGrid *theGrid,structVector *seg1, structVector *seg2);
+structTriangle* FindFirstTri(structTriangle *T, double x, double y, structGrid *theGrid, structTriangulation *theMesh, int previouTriangle,int circle);
+void splitSeg(structTriangle *Tree,structTriangulation *theTriangulation,structGrid *theGrid,structVector *seg1,structVector *seg2,int i, int a);
+void splitTri(structTriangle *Tree, structTriangulation *theTriangulation, structGrid *theGrid, double px, double py, structVector *seg1, structVector *seg2);
+void computeCircumcenter(structTriangle *Triangle,structGrid *theGrid,double *px,double *py);
+double computeAspectRatio(structTriangle *T, structGrid *theGrid);
+void readPolyFile(structGrid *theGrid, structVector *seg1, structVector *seg2,const char *inputFile, structGrid *Holes);
+int toSplit(structTriangle *T, structVector *seg1, structVector *seg2, int n, structTriangulation* theMesh, structGrid* theGrid);
+int checkAndSplitEncroachedSegments(structVector *seg1,structVector *seg2,structGrid *theGrid,structTriangulation *theMesh, structVector *segMod);
+void updateSegments(structTriangle *T, structVector *seg1, structVector *seg2, int n, structTriangulation* theMesh, structGrid* theGrid);
+void checkAndSplitSkinnyTriangle(structVector *seg1,structVector *seg2,structGrid *theGrid,structTriangulation *theTriangulation, int i);
+void removeExteriorTriangles(structTriangle *theTriangle, structTriangulation *theTriangulation, int *nTriRemoved);
+void checkOneSeg(int i, int indexCloseTri, structVector *seg1,structVector *seg2,structGrid *theGrid,structTriangulation *theMesh);
+
+
+// Quality Tree
+void addTriangle(structTriangle *theTriangle, structQualityTree *node, structTriangulation *theMesh, structQualityTree *previousL, structQualityTree *previousR);
+void update_tree(structTriangle *theTriangle, structTriangulation * theTriangulation);
+
diff --git a/initStructures.c b/initStructures.c
new file mode 100644
index 0000000..6ae3c79
--- /dev/null
+++ b/initStructures.c
@@ -0,0 +1,106 @@
+#include "header.h"
+
+
+/*
+OBJECTIFS:
+Allouer la mémoire nécessaire pour créer une structure de type triangulation.
+Une triangulation contient plusieurs parametres:
+-> nTriangles: nombre total de triangles dans la tringulation courante.
+-> Triangles: pointeurs vers les triangles de la triangulation courante.
+-> size : memoire allouee (= nombre de triangles estime)
+-> tree : pointeur vers l'arbre qui trie les triangles suivant leur qualite.
+-> begin : pointeur vers le triangle de moins bonne qualite.
+*/
+
+structTriangulation *TriangulationNew(structGrid *theGrid)
+{
+ int N = theGrid->nNode*10;
+ structTriangulation *myTriangulation = malloc(sizeof(structTriangulation));
+ myTriangulation->Triangles = malloc(sizeof(structTriangle*)*N);
+ myTriangulation->nTriangles = 0;
+ myTriangulation->size = N;
+ myTriangulation->tree = QualityTreeNew();
+ myTriangulation->begin = NULL;
+
+ return myTriangulation;
+}
+
+/*
+OBJECTIFS:
+Allouer la mémoire nécessaire pour créer une structure de type arbre.
+Un arbre contient plusieurs parametres:
+-> left : pointeur vers l'enfant de gauche (moins bonne qualite que le triangle parent)
+-> right: pointeur vers l'enfant de droite (meilleure qualite que le triangle parent)
+-> triangle : triangle associe au noeud courant de l'arbre
+-> next : pointeur vers le triangle de meilleur qualite que le triangle courant (liste chainee)
+-> quality : qualite du triangle courant
+*/
+structQualityTree * QualityTreeNew()
+{
+ structQualityTree *tree = malloc(sizeof(structQualityTree));
+ tree->left = NULL;
+ tree->right = NULL;
+ tree->next = NULL;
+ tree->quality = -1.0;
+ tree->triangle = NULL;
+ return tree;
+}
+
+/*
+OBJECTIFS:
+Allouer la mémoire nécessaire pour créer une structure de type Triangle. Un triangle
+est defini par plusieurs paremetres:
+-> Nodes : indices des 3 sommets du triangle dans la structure theGrid.
+-> Voisins : indices des 3 triangles adjacents au triangle courant.
+-> index : indice du triangle courant dans la triangulation (un indice -1
+ signifie que le triangle n'est pas dans le triangulation).
+-> flagged : boolean utile dans plusieurs fonctions:
+ >>> indique si le triangle a deja ete teste lors de la recherche de la cavite.
+ >>> indique si le triangle a deja ete supprime de la triangulation
+-> quality : qualite du triangle.
+-> area : aire du triangle.
+-> Segment : si une arete correspond a un segment, alors l'indice de ce segment est
+ stocke, sinon la valeur est -1 par defaut.
+-> toBeDeleted : boolean necessaire pour supprimer les triangles exterieurs a la fin
+ de l'algorithme.
+*/
+
+structTriangle *TriangleNew(int n1, int n2, int n3, structGrid *theGrid)
+{
+ structTriangle *myTriangle = malloc(sizeof(structTriangle));
+ myTriangle->Nodes[0] = n1;
+ myTriangle->Nodes[1] = n2;
+ myTriangle->Nodes[2] = n3;
+ myTriangle->flagged = 0;
+ myTriangle->quality = computeAspectRatio(myTriangle,theGrid);
+ myTriangle->area = computeTriangleArea(myTriangle,theGrid);
+ myTriangle->Segment[0] = -1;
+ myTriangle->Segment[1] = -1;
+ myTriangle->Segment[2] = -1;
+ myTriangle->toBeDeleted = 0;
+ return myTriangle;
+}
+
+
+/*
+OBJECTIFS:
+Allouer la memoire necessaire pour une structure de type vecteur.
+Un vecteur est defini par:
+-> nElem : nombre d'elements qu'il contient (initialement 0).
+-> size: taille de la memoire allouee.
+-> elem : vecteur d'entiers (de taille "size").
+-> lengthMinSeg: longeur du plus petit segment.
+*/
+
+
+structVector *newVector(int size)
+{
+ structVector *myVector = malloc(sizeof(structVector));
+ myVector->nElem = 0;
+ myVector->elem = malloc(sizeof(int)*size);
+ myVector->size = size;
+ myVector->lengthMinSeg = 10000;
+ return myVector;
+}
+
+
diff --git a/isInFunctions.c b/isInFunctions.c
new file mode 100644
index 0000000..602d0b6
--- /dev/null
+++ b/isInFunctions.c
@@ -0,0 +1,306 @@
+#include "header.h"
+
+/*
+OBJECTIF:
+Determiner si le point de coordonnees (x,y) se trouve dans le triangle theTriangle.
+
+INPUTS:
+-> Structure theGrid contenant les coordonnees de l'ensemble des points du maillage.
+-> Pointeur vers un triangle.
+-> Coordonnees (x,y) du point a tester.
+
+OUTPUT:
+-> Renvoie 1 si le point (x,y) se trouve a l'interieur ou sur une arete du triangle.
+-> Renvoie 0 si le point (x,y) est strictement a l'exterieur du triangle.
+*/
+
+int isInTriangle(structGrid *theGrid, structTriangle *theTriangle, double x, double y)
+{
+ double x1 = theGrid->X[theTriangle->Nodes[0]-1];
+ double y1 = theGrid->Y[theTriangle->Nodes[0]-1];
+ double x2 = theGrid->X[theTriangle->Nodes[1]-1];
+ double y2 = theGrid->Y[theTriangle->Nodes[1]-1];
+ double x3 = theGrid->X[theTriangle->Nodes[2]-1];
+ double y3 = theGrid->Y[theTriangle->Nodes[2]-1];
+
+ int b1 = sameSide(x1,y1,x2,y2,x,y,x3,y3);
+ int b2 = sameSide(x2,y2,x3,y3,x,y,x1,y1);
+ int b3 = sameSide(x3,y3,x1,y1,x,y,x2,y2);
+
+ return b1*b2*b3;
+}
+
+
+
+/*
+OBJECTIF:
+Determiner si deux points (xp,yp) et (xq,yq) sont situes du meme cote de la droite
+passant par les points (x1,y1) et (x2,y2). Une certaine tolerance est laissee pour
+traiter le cas de (xp,yp) ou (xq,yq) presque colineaires avec (x1,y1) et (x2,y2).
+
+INPUTS:
+-> Deux points (x1,y1) et (x2,y2) definissant une droite.
+-> Deux points (xp,yp) et (xq,yq) a tester.
+
+OUTPUT:
+-> Renvoie 1 si les points (xp,yp) et (xq,yq) sont situes du meme cote de la droite
+ definie par (x1,y1) et (x2,y2).
+-> Renvoie 0 sinon.
+*/
+
+int sameSide(double x1, double y1, double x2, double y2, double xp, double yp, double xq, double yq)
+{
+ double sidep = (x2-x1)*(yp-y1)-(xp-x1)*(y2-y1);
+ double sideq = (x2-x1)*(yq-y1)-(xq-x1)*(y2-y1);
+
+ if (sidep*sideq >= -0.0)
+ return 1;
+ else
+ return 0;
+}
+
+
+
+/*
+OBJECTIFS:
+Determiner si le point de coodonnees (dx,dy) se trouve a l'interieur du cercle circonscrit
+a un triangle. La caracteristique de cette fonction est de ne pas calculer directement le
+centre et le rayon du cercle mais de se baser uniquement sur les coordonnees des sommets du
+triangle.
+
+REFERENCE:
+http://www.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-predicates.ps
+
+INPUTS:
+-> Coordonnees du points a tester : (dx,dy)
+-> Pointeur vers le triangle a tester.
+-> Structure theGrid contenant les coordonnees de l'ensemble des points du maillage.
+
+OUTPUT:
+-> Renvoie 1 si le point (dx,dy) est compris a l'interieur ou sur la circonference du
+ cercle circonscrit au triangle T.
+-> Renvoie 0 sinon.
+*/
+
+int inCircle(double dx, double dy, structTriangle *T, structGrid* theGrid)
+{
+ double ax = theGrid->X[T->Nodes[0]-1];
+ double ay = theGrid->Y[T->Nodes[0]-1];
+
+ double bx = theGrid->X[T->Nodes[1]-1];
+ double by = theGrid->Y[T->Nodes[1]-1];
+
+ double cx = theGrid->X[T->Nodes[2]-1];
+ double cy = theGrid->Y[T->Nodes[2]-1];
+
+ double det = 0;
+ det -= (ax-dx)*((by-dy)*((cx-dx)*(cx-dx)+(cy-dy)*(cy-dy))-(cy-dy)*((bx-dx)*(bx-dx)+(by-dy)*(by-dy)));
+ det += (bx-dx)*((ay-dy)*((cx-dx)*(cx-dx)+(cy-dy)*(cy-dy))-(cy-dy)*((ax-dx)*(ax-dx)+(ay-dy)*(ay-dy)));
+ det -= (cx-dx)*((ay-dy)*((bx-dx)*(bx-dx)+(by-dy)*(by-dy))-(by-dy)*((ax-dx)*(ax-dx)+(ay-dy)*(ay-dy)));
+
+ if (det>0) return 1; else return 0;
+}
+
+
+/*
+OBJECTIF:
+Calculer la surface d'un triangle en se basant uniquement sur les coordonnees de ses sommets
+Cette fonction utilise la relation de HERON.
+
+INPUTS:
+-> Pointeur vers un triangle.
+-> Structure theGrid contenant les coordonnees de l'ensemble des points du maillage.
+
+OUTPUTS:
+-> Renvoie l'aire du triangle.
+*/
+
+double computeTriangleArea(structTriangle *theTriangle, structGrid *theGrid)
+{
+ // Etape 1 : recuperer les coordonnees des sommets du triangle
+
+ double x0 = theGrid->X[theTriangle->Nodes[0]-1];
+ double x1 = theGrid->X[theTriangle->Nodes[1]-1];
+ double x2 = theGrid->X[theTriangle->Nodes[2]-1];
+
+ double y0 = theGrid->Y[theTriangle->Nodes[0]-1];
+ double y1 = theGrid->Y[theTriangle->Nodes[1]-1];
+ double y2 = theGrid->Y[theTriangle->Nodes[2]-1];
+
+ // Etape 2 : calculer la longueur des aretes du triangle
+
+ double a = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
+ double b = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0));
+ double c = sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1));
+
+ // Etape 3 : calculer le demi-perimetre et la surface via la formule de HERON
+
+ double s = 0.5*(a+b+c);
+ double S = sqrt(s*(s-a)*(s-b)*(s-c));
+
+ return S;
+}
+
+
+/*
+OBJECTIF:
+Determiner si le point de coordonnees (cx,cy) se trouve dans dans le cercle diametral du
+segment defini par les sommets (ax,ay) et (bx,by).
+
+INPUTS:
+-> Coordonnees des deux sommets du segment
+-> Coordonnees du point a tester
+
+OUTPUTS:
+-> Renvoie 1 si le point (cx,cy) se trouve dans le cercle diametral du segment defini par
+ les sommets (ax,ay) et (bx,by). Renvoie 0 sinon.
+*/
+
+int inDiametralCircle(double ax, double ay, double bx, double by, double cx, double cy)
+{
+ double centerX = 0.5*(ax+bx);
+ double centerY = 0.5*(ay+by);
+ double Rsq = sqrt((centerX-ax)*(centerX-ax) + (centerY-ay)*(centerY-ay));
+ double distsq = sqrt((centerX-cx)*(centerX-cx) + (centerY-cy)*(centerY-cy));
+
+ if (distsq<=Rsq)
+ return 1;
+ else
+ return 0;
+}
+
+
+/*
+OBJECTIF:
+Determiner si le point de coordonnees (x,y) se trouve dans dans le cercle diametral du
+segment de numero "ind".
+
+INPUTS:
+-> Vecteurs definissant les extremites de l ensemble des segments
+-> Structure theGrid contenant les coordonnees de l'ensemble des points du maillage.
+-> Coordonnees du point a tester
+-> Numero du segment a tester.
+
+OUTPUTS:
+-> Renvoie 1 si le point (x,y) se trouve dans le cercle diametral du segment teste.
+ Renvoie 0 sinon.
+*/
+int checkSeg(structVector *seg1,structVector *seg2,structGrid *theGrid, double x, double y, int ind)
+{
+ int n1 = seg1->elem[ind]-1;
+ int n2 = seg2->elem[ind]-1;
+ double ax = theGrid->X[n1];
+ double ay = theGrid->Y[n1];
+ double bx = theGrid->X[n2];
+ double by = theGrid->Y[n2];
+
+ return inDiametralCircle(ax, ay, bx, by, x, y);
+}
+
+
+/*
+OBJECTIF:
+Calculer l aspect ratio d un triangle (rapport entre le rayon du cercle circonscrit au
+triangle et la plus petite arete de ce triangle)
+
+INPUTS:
+-> Pointeur vers un triangle
+-> Structure theGrid contenant les coordonnees de l'ensemble des points du maillage.
+
+OUTPUTS:
+-> Renvoie l'aspect ratio du triangle.
+*/
+double computeAspectRatio(structTriangle *T, structGrid *theGrid)
+{
+
+ double x0 = theGrid->X[T->Nodes[0]-1];
+ double x1 = theGrid->X[T->Nodes[1]-1];
+ double x2 = theGrid->X[T->Nodes[2]-1];
+
+ double y0 = theGrid->Y[T->Nodes[0]-1];
+ double y1 = theGrid->Y[T->Nodes[1]-1];
+ double y2 = theGrid->Y[T->Nodes[2]-1];
+
+ double a = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
+ double b = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0));
+ double c = sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1));
+
+ double shortestEdge = Min(a,Min(b,c));
+
+ double s = 0.5*(a+b+c);
+ double S = sqrt(s*(s-a)*(s-b)*(s-c)); // Heron's formula
+ double radius = (a*b*c)/(4.0*S);
+
+ extern double qualityRequest;
+ extern double areaRequest;
+
+ double area = computeTriangleArea(T,theGrid);
+
+ return Max(radius/shortestEdge/qualityRequest,area/areaRequest);
+}
+
+
+/*
+OBJECTIF:
+Calculer les coordonnees du centre du cercle circonscrit a un triangle.
+
+INPUTS:
+-> Pointeur vers un triangle
+-> Structure theGrid contenant les coordonnees de l'ensemble des points du maillage.
+-> Pointeurs vers les coordonnees du centre du cercle.
+
+OUTPUTS:
+-> Les coordonnees du centre du cercle sont mises a jour.
+*/
+void computeCircumcenter(structTriangle *Triangle,structGrid *theGrid,double *px,double *py)
+{
+ double x0 = theGrid->X[Triangle->Nodes[0]-1];
+ double x1 = theGrid->X[Triangle->Nodes[1]-1];
+ double x2 = theGrid->X[Triangle->Nodes[2]-1];
+
+ double y0 = theGrid->Y[Triangle->Nodes[0]-1];
+ double y1 = theGrid->Y[Triangle->Nodes[1]-1];
+ double y2 = theGrid->Y[Triangle->Nodes[2]-1];
+
+ // Compute the center
+ double p1,p2,m1,m2;
+ if(fabs(y2-y1)>=fabs(y1-y0) && fabs(y2-y0)>=fabs(y1-y0))
+ {
+ m1 = - (x2-x1)/(y2-y1);
+ m2 = - (x2-x0)/(y2-y0);
+ p1 = (y2+y1)*0.5 - m1*(x2+x1)*0.5;
+ p2 = (y2+y0)*0.5 - m2*(x2+x0)*0.5;
+ }
+ else if(fabs(y1-y0)>=fabs(y2-y1) && fabs(y2-y0)>=fabs(y2-y1))
+ {
+ m1 = - (x1-x0)/(y1-y0);
+ m2 = - (x2-x0)/(y2-y0);
+ p1 = (y1+y0)*0.5 - m1*(x1+x0)*0.5;
+ p2 = (y2+y0)*0.5 - m2*(x2+x0)*0.5;
+
+ }
+ else if(fabs(y1-y0)>=fabs(y2-y0) && fabs(y2-y1)>=fabs(y2-y0))
+ {
+ m1 = - (x1-x0)/(y1-y0);
+ m2 = - (x2-x1)/(y2-y1);
+ p1 = (y1+y0)*0.5 - m1*(x1+x0)*0.5;
+ p2 = (y2+y1)*0.5 - m2*(x2+x1)*0.5;
+
+ }
+ else
+ {
+ printf("ATTENTION DIVISION PAR 0 - CENTRE CERCLE !!! \n");
+ }
+
+ double deltaM = m1-m2;
+
+ *px = (p2-p1)/deltaM;
+ *py = m1*(*px)+p1;
+
+ if(!inCircle(*px,*py, Triangle, theGrid))
+ {
+ printf("Nodes: y0 %f y1 %f y2 %f \n",y0,y1,y2);
+ printf("Error computing the circumcenter\n");
+ exit(0);
+ }
+}
diff --git a/main.c b/main.c
new file mode 100644
index 0000000..13e7ade
--- /dev/null
+++ b/main.c
@@ -0,0 +1,175 @@
+#include "header.h"
+
+double qualityRequest;
+double areaRequest;
+
+int main(int argc, char *argv[])
+{
+
+ // Lancement du programme : l'utilisateur peut demander un message d'aide (tapper -help).
+ // Par defaut le programme est lance sur le fichier superior.poly. L'utilisateur peut
+ // egalement choisir un autre fichier contenant un nuage de points.
+
+ char *fileName; double angle;
+ char file_directory[100] = "";
+ if (argc == 1)
+ {
+ fileName = "../superior.poly";
+ angle = 20.7;
+ areaRequest = 1e14;
+ }
+ else if (argc == 2)
+ {
+ if (!strcmp(argv[1],"-help"))
+ {
+ printManPage("../help.txt");
+ exit(0);
+ }
+ else
+ {
+ fileName = strcat(file_directory,argv[1]);
+ angle = 20.7;
+ areaRequest = 1e14;
+ }
+ }
+ else if (argc == 3)
+ {
+ fileName = strcat(file_directory,argv[1]);
+ angle = atof(argv[2]);
+ areaRequest = 1e14;
+ }
+ else if (argc == 4)
+ {
+ fileName = strcat(file_directory,argv[1]);
+ angle = atof(argv[2]);
+ areaRequest = atof(argv[3])*100;
+ }
+ else
+ {
+ printf("Error number arguments");
+ }
+
+ clock_t start;
+
+ // Etape 1 : Demander a l utilisateur les parametres voulus : angle minimal et
+ // arbre pour classer les triangles suivant leur qualite
+
+ double minAngle = angle*3.1415/180.0;
+ qualityRequest = 1.0/(2.0*sin(minAngle));
+
+
+
+
+ // Etape 2 : Initialisation des structures et lecture du fichier de donnees
+
+ structGrid *theGrid = malloc(sizeof(structGrid));
+ structGrid *Holes = malloc(sizeof(structGrid));
+ structVector *seg1 = newVector(10000);
+ structVector *seg2 = newVector(10000);
+ readPolyFile(theGrid,seg1, seg2,fileName,Holes);
+ structTriangulation *theTriangulation = TriangulationNew(theGrid);
+
+ if(theTriangulation->tree->left != NULL || theTriangulation->tree->right != NULL)
+ {
+ printf("Error! \n"); exit(0);
+ }
+
+
+ // Etape 3: Ajout d un grand carre autour de l ensemble de points donnes
+
+ addExteriorNodes(theGrid,seg1,seg2,theTriangulation);
+
+
+
+
+ // Etape 4 : Division des segments qui contiennent un sommet dans leur cercle diametral
+
+ start = clock();
+ printf("\n\nInput data:");
+ printf("\n -> Number of points: %d ",theGrid->nNode);
+ printf("\n -> Number of segments: %d \n",seg1->nElem);
+ checkAndSplitEncroachedSegments(seg1,seg2,theGrid,theTriangulation,NULL);
+
+ //printf("\n\nEncroached segments divided:");
+ //printf("\n -> Number of points: %d ",theGrid->nNode);
+ //printf("\n -> Number of triangles: %d",theTriangulation->nTriangles);
+ //printf("\n -> Time: %f [sec]\n\n",(double)(clock() - start)/CLOCKS_PER_SEC);
+
+
+ // Etape 5 : Division des mauvais triangles (en mettant tour a tour les encroached segments a jour)
+ start = clock();
+ while(1)
+ {
+ structTriangle *T = theTriangulation->begin->triangle;
+
+ // Si l'arbre contient un triangle de mauvaise qualite
+ if(T->quality>1)
+ {
+ // Si le triangle n'existe plus
+ if(T->flagged)
+ {
+ // Mettre une tres mauvaise qualite afin de ne plus rajouter de triangle a gauche du triangle supprime
+ theTriangulation->begin->quality = 1e14; free(T); theTriangulation->begin->triangle = NULL;
+ theTriangulation->begin = theTriangulation->begin->next;
+ continue;
+ }
+
+ // Si le triangle existe, il doit etre divise (et on modifie prealablement les segments si ils en sont la cause)
+ canBeSplit(theTriangulation->Triangles[T->index],seg1, seg2, theGrid, minAngle,theTriangulation);
+ checkAndSplitSkinnyTriangle(seg1,seg2,theGrid,theTriangulation,T->index);
+ }
+ else
+ {
+ break;
+ }
+ }
+
+
+ double myQuality = theTriangulation->begin->quality;
+
+ printf("\n\nRuppert's algorithm finished:");
+ printf("\n -> Number of points: %d ",theGrid->nNode);
+ printf("\n -> Number of triangles: %d",theTriangulation->nTriangles);
+ printf("\n -> Minimum angle > %f [deg]",asin(1.0/(2.0*myQuality*qualityRequest))*180.0/3.1415);
+ printf("\n -> Time: %f [sec] \n",(double)(clock() - start)/CLOCKS_PER_SEC);
+ showIt(theGrid, theTriangulation, seg1, seg2, 1, 0,angle);
+
+
+ // Etape 6 : Enlever les triangles exterieurs
+
+ start = clock();
+ int nTriRemoved = 0;
+ double x = theGrid->X[theGrid->nInputVertex-4]+1e-12;
+ double y = theGrid->Y[theGrid->nInputVertex-4]+1e-12;
+ structTriangle* firstTRI = FindFirstTri(theTriangulation->Triangles[0], x, y, theGrid, theTriangulation, -1, 0);
+ removeExteriorTriangles(firstTRI,theTriangulation,&nTriRemoved);
+ int k=0;
+ for(k=0;k<Holes->nNode;k++)
+ {
+ firstTRI = FindFirstTri(theTriangulation->Triangles[0], Holes->X[k], Holes->Y[k], theGrid, theTriangulation, -1, 0);
+ removeExteriorTriangles(firstTRI,theTriangulation,&nTriRemoved);
+ }
+
+ //printf("\n\nRemoving exterior triangles:");
+ //printf("\n -> Number of triangles removed: %d",nTriRemoved);
+ //printf("\n -> Time: %f [sec]",(double)(clock() - start)/CLOCKS_PER_SEC);
+ //printf("\n\n");
+
+
+ showIt(theGrid, theTriangulation, seg1, seg2, 1, 0,angle);
+
+
+ // Etape 7 : Liberation de la mémoire
+
+ freeQualityTree(theTriangulation->tree);
+ freeVector(seg1);
+ freeVector(seg2);
+ freeTriangulation(theTriangulation);
+ freeGrid(theGrid);
+ freeGrid(Holes);
+}
+
+
+
+
+
diff --git a/makefile.linux b/makefile.linux
new file mode 100644
index 0000000..78fec3c
--- /dev/null
+++ b/makefile.linux
@@ -0,0 +1,7 @@
+COMPILER = gcc -pg
+OBJECTS = *.c
+OUTPUT = ruppertL
+LIBS = $(SUBLIBS) -lglfw -lXext -lX11 -lm -lpthread -lGL -lGLU
+
+$(OUTPUT) : $(OBJECTS)
+ $(COMPILER) $(OBJECTS) -o $(OUTPUT) $(LIBS)
diff --git a/printManPage.c b/printManPage.c
new file mode 100644
index 0000000..040e861
--- /dev/null
+++ b/printManPage.c
@@ -0,0 +1,22 @@
+#include "header.h"
+
+void printManPage(const char *inputFile)
+{
+ FILE* fichier = fopen(inputFile,"r");
+ char chaine[100] = "";
+
+ if (fichier != NULL)
+ {
+ while (fgets(chaine, 100, fichier) != NULL) // On lit le fichier tant qu'on ne reçoit pas d'erreur (NULL)
+ {
+ printf("%s", chaine); // On affiche la chaîne qu'on vient de lire
+ }
+
+ fclose(fichier);
+ }
+ else
+ {
+ printf("Erreur à la lecture du help \n");
+ }
+
+}
diff --git a/readInputFile.c b/readInputFile.c
new file mode 100644
index 0000000..c66468b
--- /dev/null
+++ b/readInputFile.c
@@ -0,0 +1,150 @@
+#include "header.h"
+
+void readPolyFile(structGrid *theGrid, structVector *seg1, structVector *seg2,const char *inputFile, structGrid *Holes)
+{
+ FILE* file = fopen(inputFile,"r");
+ int i,nVertex,dim,nAttributes,nBoundaryMarkers,nSeg, nBoundaryMarkers2,nHoles;
+ double number,x,y,attribute,boundary,ext1,ext2;
+
+ if(file!=NULL)
+ {
+ // Etape 1 : lecture des coordonnees des sommets
+
+ fscanf(file,"%d %d %d %d \n",&nVertex,&dim,&nAttributes,&nBoundaryMarkers);
+ theGrid->nNode = nVertex;
+ theGrid->size = theGrid->nNode*10;
+
+ theGrid->X = malloc(sizeof(double)*theGrid->size);
+ theGrid->Y = malloc(sizeof(double)*theGrid->size);
+ theGrid->Z = malloc(sizeof(double)*theGrid->size);
+
+ for (i = 0; i < theGrid->nNode; ++i)
+ {
+ if(nBoundaryMarkers==0)
+ {
+ if(nAttributes==0)
+ fscanf(file,"%le %le %le \n",&number,&x, &y);
+ else
+ fscanf(file,"%le %le %le %le\n",&number,&x, &y, &attribute);
+ }
+ else
+ {
+ if(nAttributes==0)
+ fscanf(file,"%le %le %le %le \n",&number,&x, &y,&boundary);
+ else
+ fscanf(file,"%le %le %le %le %le \n",&number,&x, &y, &attribute,&boundary);
+ }
+
+ double r1 = rand() / (double)RAND_MAX;
+ double r2 = rand() / (double)RAND_MAX;
+
+ theGrid->X[i] = x+r1*1e-11;
+ theGrid->Y[i] = y+r2*1e-11;
+ theGrid->Z[i] = 0.0;
+ }
+
+
+ // Etape 2 : lecture des extremites des segments.
+
+ fscanf(file,"%d %d \n",&nSeg,&nBoundaryMarkers2);
+ seg1->nElem = nSeg;
+ seg2->nElem = nSeg;
+
+ if(seg1->nElem>=seg1->size)
+ {
+ printf("\n\n----------------------------------------------------- \n");
+ printf("Realloc of seg1->elem and seg2->elem \n");
+ printf("----------------------------------------------------- \n");
+ int N = seg1->nElem*100;
+ seg1->elem = realloc(seg1->elem,N*sizeof(int));
+ seg2->elem = realloc(seg2->elem,N*sizeof(int));
+ seg1->size = seg1->size*N;
+ seg2->size = seg1->size*N;
+ }
+
+ for (i = 0; i < nSeg; ++i)
+ {
+ if(nBoundaryMarkers2==0)
+ fscanf(file,"%le %le %le \n",&number,&ext1, &ext2);
+ else
+ fscanf(file,"%le %le %le %le \n",&number,&ext1, &ext2,&boundary);
+
+ seg1->elem[i] = ext1;
+ seg2->elem[i] = ext2;
+ }
+
+ // Etape 3 : lecture des "trous"
+
+ fscanf(file,"%d \n",&nHoles);
+ Holes->nNode = nHoles;
+ Holes->size = nHoles;
+ Holes->X = malloc(sizeof(double)*Holes->size);
+ Holes->Y = malloc(sizeof(double)*Holes->size);
+ Holes->Z = malloc(sizeof(double)*Holes->size);
+ for(i=0; i<nHoles; i++)
+ {
+ fscanf(file,"%le %le %le \n",&number,&x, &y);
+ Holes->X[i] = x;
+ Holes->Y[i] = y;
+ }
+
+ // Etape 4 : le nuage de points est centre en (0,0) et les donnees sont
+ // mise a l'echelle pour faciliter le zoom.
+
+ double* MyX = theGrid->X;
+ double* MyY = theGrid->Y;
+
+ double minX = min(MyX,theGrid->nNode);
+ double maxX = max(MyX,theGrid->nNode);
+ double minY = min(MyY,theGrid->nNode);
+ double maxY = max(MyY,theGrid->nNode);
+
+ double shiftX = -(maxX+minX)/2.0;
+ double shiftY = -(maxY+minY)/2.0;
+
+ int l;
+ for(l=0;l<theGrid->nNode;l++)
+ {
+ theGrid->X[l] += shiftX;
+ theGrid->Y[l] += shiftY;
+ }
+ for(l=0;l<Holes->nNode;l++)
+ {
+ Holes->X[l] += shiftX;
+ Holes->Y[l] += shiftY;
+ }
+ for(l=0;l<theGrid->nNode;l++)
+ {
+ theGrid->X[l] *= 500/(maxX-minX+maxY-minY);
+ theGrid->Y[l] *= 500/(maxX-minX+maxY-minY);
+ }
+ for(l=0;l<Holes->nNode;l++)
+ {
+ Holes->X[l] *= 500/(maxX-minX+maxY-minY);
+ Holes->Y[l] *= 500/(maxX-minX+maxY-minY);
+ }
+
+ // Etape 5 : calcul de la longueur du plus petit segment.
+
+ double lengthMinSeg = 100000;
+ for (i = 0; i < nSeg; ++i)
+ {
+ double ax = theGrid->X[seg1->elem[i]-1];
+ double bx = theGrid->X[seg2->elem[i]-1];
+ double ay = theGrid->Y[seg1->elem[i]-1];
+ double by = theGrid->Y[seg2->elem[i]-1];
+ double length = sqrt((ax-bx)*(ax-bx) + (ay-by)*(ay-by));
+ lengthMinSeg = Min(lengthMinSeg,length);
+ }
+ seg1->lengthMinSeg = lengthMinSeg;
+ seg2->lengthMinSeg = lengthMinSeg;
+
+ }
+ else
+ {
+ printf("ERROR in function readInputFile : file = NULL \n");
+ exit(0);
+ }
+ fclose(file);
+}
+
diff --git a/removeExteriorTriangles.c b/removeExteriorTriangles.c
new file mode 100644
index 0000000..63a04dd
--- /dev/null
+++ b/removeExteriorTriangles.c
@@ -0,0 +1,39 @@
+#include "header.h"
+
+/*
+OBJECTIF:
+Une fois que la triangulation est terminee, supprimer les triangles situes a
+l'exterieur de la forme a mailler. Technique recursive.
+
+INPUTS:
+-> theTriangle : pointeur vers un des triangles a supprimer.
+-> theTriangulation : pointeur vers l'ensemble des triangles de la triangulation
+ courante.
+-> nTriRemoved : initialment 0
+
+OUTPUTS:
+-> Tous les triangles situes a l'exterieur ont ete marques "toBeDeleted = 1"
+*/
+void removeExteriorTriangles(structTriangle *theTriangle, structTriangulation *theTriangulation, int *nTriRemoved)
+{
+ // Etape 1 : le triangle courant est marque "toBeDeleted"
+
+ theTriangle->toBeDeleted = 1;
+ *nTriRemoved = *nTriRemoved + 1;
+
+ // Etape 2 : parcourir les voisins du triangle courant. Un voisin ne sera
+ // marque "toBeDeleted" que si ce voisin n'a pas encore ete teste et si
+ // l'arete entre le triangle et ce voisin n'est pas un segment.
+
+ int i = 0;
+ for(i=0;i<3;i++)
+ {
+ if(theTriangle->Segment[i]==-1)
+ {
+ if(theTriangulation->Triangles[theTriangle->Voisins[i]]->toBeDeleted==0)
+ {
+ removeExteriorTriangles(theTriangulation->Triangles[theTriangle->Voisins[i]],theTriangulation,nTriRemoved);
+ }
+ }
+ }
+}
diff --git a/superior.poly b/superior.poly
new file mode 100644
index 0000000..0b96356
--- /dev/null
+++ b/superior.poly
@@ -0,0 +1,1045 @@
+518 2 0 0
+1 0.1 -0.629607
+2 0.102293 -0.635353
+3 0.107467 -0.638807
+4 0.11092 -0.643967
+5 0.12184 -0.643967
+6 0.128733 -0.643393
+7 0.135053 -0.641673
+8 0.141373 -0.63938
+9 0.147707 -0.637073
+10 0.154027 -0.63478
+11 0.164373 -0.633633
+12 0.171267 -0.633047
+13 0.177013 -0.630753
+14 0.18104 -0.626153
+15 0.18908 -0.62386
+16 0.194253 -0.626727
+17 0.200573 -0.62558
+18 0.205747 -0.622713
+19 0.210347 -0.618687
+20 0.214373 -0.614087
+21 0.221267 -0.611207
+22 0.227013 -0.608913
+23 0.23276 -0.61006
+24 0.236787 -0.61466
+25 0.239653 -0.619833
+26 0.237933 -0.62558
+27 0.235053 -0.631327
+28 0.230467 -0.635353
+29 0.225867 -0.639953
+30 0.22816 -0.6457
+31 0.225293 -0.651447
+32 0.221267 -0.656047
+33 0.218387 -0.66122
+34 0.222987 -0.665233
+35 0.22816 -0.661793
+36 0.233333 -0.658913
+37 0.23908 -0.655473
+38 0.244827 -0.65374
+39 0.249427 -0.657193
+40 0.254027 -0.66122
+41 0.258613 -0.66466
+42 0.264373 -0.665233
+43 0.27012 -0.668113
+44 0.275867 -0.668113
+45 0.282187 -0.665233
+46 0.28908 -0.662367
+47 0.294827 -0.66006
+48 0.302293 -0.657193
+49 0.30804 -0.656047
+50 0.3138 -0.654887
+51 0.318387 -0.650873
+52 0.322987 -0.646273
+53 0.327587 -0.642247
+54 0.332187 -0.637647
+55 0.340227 -0.63478
+56 0.35 -0.633633
+57 0.355747 -0.635927
+58 0.361493 -0.63478
+59 0.36896 -0.633047
+60 0.375867 -0.631327
+61 0.38104 -0.62846
+62 0.386787 -0.625007
+63 0.392533 -0.621553
+64 0.397707 -0.618113
+65 0.402293 -0.61466
+66 0.407467 -0.610633
+67 0.414373 -0.607767
+68 0.42012 -0.606047
+69 0.42644 -0.60374
+70 0.429307 -0.597993
+71 0.43276 -0.593393
+72 0.437933 -0.589953
+73 0.442533 -0.5865
+74 0.447133 -0.582473
+75 0.45172 -0.577887
+76 0.45632 -0.574433
+77 0.461493 -0.57098
+78 0.46724 -0.56754
+79 0.471267 -0.56294
+80 0.475867 -0.557767
+81 0.480467 -0.553167
+82 0.486787 -0.549713
+83 0.49196 -0.546847
+84 0.49828 -0.543967
+85 0.5046 -0.5411
+86 0.510347 -0.538793
+87 0.516093 -0.5365
+88 0.52356 -0.535353
+89 0.53104 -0.5365
+90 0.540227 -0.537647
+91 0.546547 -0.539953
+92 0.552867 -0.542247
+93 0.541373 -0.546847
+94 0.535627 -0.546847
+95 0.529307 -0.54914
+96 0.52816 -0.554887
+97 0.524707 -0.56006
+98 0.519533 -0.563513
+99 0.514373 -0.566967
+100 0.509773 -0.57098
+101 0.505173 -0.575007
+102 0.500573 -0.579033
+103 0.49712 -0.583633
+104 0.49368 -0.588793
+105 0.490227 -0.593967
+106 0.486787 -0.59914
+107 0.483333 -0.60374
+108 0.47988 -0.609487
+109 0.477013 -0.615807
+110 0.475293 -0.621553
+111 0.47356 -0.6273
+112 0.472987 -0.633047
+113 0.472413 -0.638807
+114 0.472987 -0.64454
+115 0.477587 -0.640527
+116 0.481613 -0.635927
+117 0.486213 -0.6319
+118 0.4908 -0.62846
+119 0.4954 -0.624433
+120 0.498853 -0.61926
+121 0.503453 -0.615807
+122 0.5092 -0.61294
+123 0.505747 -0.618113
+124 0.50172 -0.623287
+125 0.49828 -0.627887
+126 0.4954 -0.633633
+127 0.5 -0.62846
+128 0.504027 -0.623287
+129 0.509773 -0.621553
+130 0.51896 -0.62098
+131 0.52988 -0.622127
+132 0.538507 -0.623287
+133 0.54368 -0.626153
+134 0.549427 -0.629033
+135 0.555747 -0.631327
+136 0.562067 -0.633633
+137 0.56724 -0.6365
+138 0.567813 -0.642247
+139 0.57012 -0.647993
+140 0.572987 -0.65374
+141 0.575867 -0.66006
+142 0.57816 -0.66638
+143 0.582187 -0.67098
+144 0.586213 -0.67558
+145 0.587933 -0.681327
+146 0.605747 -0.682473
+147 0.618387 -0.6819
+148 0.62472 -0.680753
+149 0.630453 -0.6819
+150 0.633907 -0.6865
+151 0.638507 -0.689953
+152 0.647133 -0.691673
+153 0.64828 -0.685927
+154 0.6546 -0.685353
+155 0.6592 -0.689953
+156 0.66264 -0.69454
+157 0.668387 -0.691673
+158 0.67184 -0.6865
+159 0.675293 -0.6819
+160 0.678733 -0.6773
+161 0.685627 -0.673287
+162 0.691373 -0.669833
+163 0.695973 -0.66638
+164 0.70172 -0.66294
+165 0.70632 -0.659487
+166 0.711493 -0.656047
+167 0.71724 -0.65834
+168 0.724133 -0.65662
+169 0.732187 -0.654887
+170 0.737933 -0.65374
+171 0.74368 -0.652593
+172 0.75 -0.653167
+173 0.75804 -0.65374
+174 0.78448 -0.653167
+175 0.789653 -0.6503
+176 0.794827 -0.64742
+177 0.800573 -0.64454
+178 0.805747 -0.641673
+179 0.813213 -0.63938
+180 0.81896 -0.63938
+181 0.824707 -0.639953
+182 0.832187 -0.63822
+183 0.838507 -0.637647
+184 0.833907 -0.642247
+185 0.83104 -0.64742
+186 0.83104 -0.653167
+187 0.83104 -0.658913
+188 0.831613 -0.66466
+189 0.83104 -0.670407
+190 0.82988 -0.676727
+191 0.844253 -0.68018
+192 0.850573 -0.682473
+193 0.858627 -0.682473
+194 0.864373 -0.680753
+195 0.87012 -0.679033
+196 0.87644 -0.680753
+197 0.878733 -0.6865
+198 0.886787 -0.6865
+199 0.890227 -0.6819
+200 0.885053 -0.67846
+201 0.881613 -0.673287
+202 0.87816 -0.668687
+203 0.87988 -0.66294
+204 0.883333 -0.65834
+205 0.886787 -0.65374
+206 0.8908 -0.64914
+207 0.891373 -0.643393
+208 0.888507 -0.637073
+209 0.8862 -0.64282
+210 0.881613 -0.63938
+211 0.881613 -0.633047
+212 0.881613 -0.6273
+213 0.88736 -0.62558
+214 0.893107 -0.62558
+215 0.898853 -0.622127
+216 0.9 -0.61638
+217 0.8954 -0.61294
+218 0.888507 -0.608913
+219 0.88276 -0.60834
+220 0.877013 -0.607767
+221 0.87184 -0.611207
+222 0.866667 -0.60834
+223 0.862067 -0.604313
+224 0.855173 -0.604313
+225 0.855747 -0.597993
+226 0.85804 -0.591673
+227 0.859773 -0.585927
+228 0.860347 -0.58018
+229 0.86092 -0.574433
+230 0.864373 -0.569833
+231 0.86724 -0.56466
+232 0.870693 -0.558913
+233 0.870693 -0.553167
+234 0.86724 -0.547993
+235 0.86092 -0.545127
+236 0.859773 -0.53938
+237 0.855173 -0.53478
+238 0.85 -0.5319
+239 0.844827 -0.529033
+240 0.839653 -0.526153
+241 0.833907 -0.522713
+242 0.82988 -0.516967
+243 0.82644 -0.512367
+244 0.82816 -0.506047
+245 0.83104 -0.499713
+246 0.833907 -0.493967
+247 0.836787 -0.487647
+248 0.83908 -0.481327
+249 0.835053 -0.476727
+250 0.837933 -0.470407
+251 0.8408 -0.465233
+252 0.842533 -0.458913
+253 0.82184 -0.458913
+254 0.80632 -0.459487
+255 0.795973 -0.46006
+256 0.790227 -0.461793
+257 0.780467 -0.46294
+258 0.762067 -0.463513
+259 0.750573 -0.46294
+260 0.7454 -0.46006
+261 0.740227 -0.45662
+262 0.735053 -0.453167
+263 0.73104 -0.447993
+264 0.727587 -0.443393
+265 0.72356 -0.43822
+266 0.720693 -0.4319
+267 0.718387 -0.42558
+268 0.716093 -0.41926
+269 0.713787 -0.413513
+270 0.711493 -0.407767
+271 0.7092 -0.401447
+272 0.70632 -0.3957
+273 0.7046 -0.389953
+274 0.70172 -0.383633
+275 0.699427 -0.377887
+276 0.697133 -0.372127
+277 0.694827 -0.365807
+278 0.69196 -0.359487
+279 0.68908 -0.353167
+280 0.6862 -0.347993
+281 0.683333 -0.34282
+282 0.67644 -0.340427
+283 0.671267 -0.343393
+284 0.658627 -0.343967
+285 0.6592 -0.337647
+286 0.654027 -0.333047
+287 0.64828 -0.335353
+288 0.64368 -0.33938
+289 0.633907 -0.342247
+290 0.62816 -0.33938
+291 0.62644 -0.333047
+292 0.621267 -0.3365
+293 0.612067 -0.339953
+294 0.605173 -0.340527
+295 0.597707 -0.339953
+296 0.592533 -0.3365
+297 0.58736 -0.333633
+298 0.581613 -0.33018
+299 0.574707 -0.3273
+300 0.56896 -0.325007
+301 0.563213 -0.322713
+302 0.557467 -0.320407
+303 0.551147 -0.318113
+304 0.5454 -0.315807
+305 0.539653 -0.31466
+306 0.53276 -0.313513
+307 0.527587 -0.309487
+308 0.522413 -0.30546
+309 0.516093 -0.307767
+310 0.510347 -0.309487
+311 0.5046 -0.31006
+312 0.498853 -0.307193
+313 0.498853 -0.313513
+314 0.498853 -0.319833
+315 0.5 -0.326153
+316 0.502293 -0.3319
+317 0.50804 -0.334207
+318 0.5138 -0.337073
+319 0.5138 -0.34282
+320 0.5138 -0.348567
+321 0.513213 -0.354313
+322 0.509773 -0.358913
+323 0.505747 -0.36466
+324 0.502293 -0.36926
+325 0.498853 -0.37386
+326 0.493107 -0.373287
+327 0.488507 -0.3773
+328 0.48276 -0.38018
+329 0.47816 -0.384207
+330 0.475293 -0.38938
+331 0.469533 -0.3911
+332 0.46552 -0.385927
+333 0.468387 -0.379607
+334 0.472413 -0.375007
+335 0.477013 -0.370407
+336 0.48104 -0.365807
+337 0.485627 -0.36122
+338 0.489653 -0.35662
+339 0.492533 -0.3503
+340 0.494827 -0.343967
+341 0.49368 -0.337647
+342 0.488507 -0.333047
+343 0.480467 -0.331327
+344 0.475293 -0.335353
+345 0.472413 -0.341673
+346 0.470693 -0.347993
+347 0.470693 -0.35374
+348 0.46896 -0.359487
+349 0.4638 -0.362367
+350 0.46092 -0.36754
+351 0.4592 -0.37386
+352 0.457467 -0.379607
+353 0.455747 -0.385927
+354 0.454027 -0.391673
+355 0.452293 -0.39742
+356 0.45 -0.40374
+357 0.444253 -0.405473
+358 0.43736 -0.407193
+359 0.433907 -0.40202
+360 0.440227 -0.39742
+361 0.441373 -0.391673
+362 0.44196 -0.385927
+363 0.445973 -0.381327
+364 0.45 -0.376727
+365 0.450573 -0.37098
+366 0.44368 -0.370407
+367 0.437933 -0.372713
+368 0.431613 -0.375007
+369 0.425867 -0.376727
+370 0.417813 -0.379033
+371 0.412067 -0.380753
+372 0.406893 -0.385353
+373 0.40288 -0.389953
+374 0.402293 -0.396273
+375 0.40288 -0.40202
+376 0.403453 -0.40834
+377 0.401147 -0.41466
+378 0.39828 -0.42098
+379 0.395973 -0.4273
+380 0.39368 -0.433047
+381 0.389653 -0.437647
+382 0.385053 -0.4411
+383 0.380467 -0.444553
+384 0.375293 -0.44742
+385 0.37012 -0.451447
+386 0.36552 -0.454887
+387 0.360347 -0.457767
+388 0.355747 -0.46122
+389 0.351147 -0.46466
+390 0.345973 -0.468113
+391 0.3408 -0.47098
+392 0.33448 -0.47386
+393 0.32816 -0.476727
+394 0.322987 -0.479607
+395 0.316093 -0.482473
+396 0.310347 -0.484207
+397 0.3046 -0.485927
+398 0.298853 -0.487647
+399 0.292533 -0.48938
+400 0.286787 -0.4911
+401 0.28104 -0.493967
+402 0.274707 -0.496847
+403 0.268387 -0.499713
+404 0.26264 -0.502593
+405 0.25632 -0.50546
+406 0.25 -0.50834
+407 0.244827 -0.511207
+408 0.23908 -0.51466
+409 0.23448 -0.51926
+410 0.229307 -0.523287
+411 0.224707 -0.5273
+412 0.219533 -0.5319
+413 0.214947 -0.535927
+414 0.210347 -0.539953
+415 0.205173 -0.543967
+416 0.2 -0.548567
+417 0.1954 -0.553167
+418 0.1908 -0.55662
+419 0.185627 -0.561207
+420 0.181613 -0.565807
+421 0.177013 -0.570407
+422 0.172413 -0.574433
+423 0.16724 -0.579033
+424 0.162067 -0.583633
+425 0.156893 -0.58822
+426 0.15172 -0.592247
+427 0.146547 -0.596847
+428 0.14196 -0.6003
+429 0.136787 -0.604887
+430 0.132187 -0.60834
+431 0.125867 -0.611207
+432 0.120693 -0.614087
+433 0.114373 -0.616967
+434 0.1092 -0.619833
+435 0.105173 -0.624433
+436 0.237933 -0.641673
+437 0.23448 -0.637073
+438 0.238507 -0.6319
+439 0.24368 -0.629033
+440 0.24828 -0.62558
+441 0.25288 -0.621553
+442 0.25804 -0.625007
+443 0.25288 -0.627887
+444 0.246547 -0.629607
+445 0.246547 -0.635353
+446 0.2408 -0.6365
+447 0.251147 -0.616967
+448 0.2546 -0.611793
+449 0.262067 -0.608913
+450 0.26092 -0.61466
+451 0.255173 -0.61638
+452 0.27356 -0.602593
+453 0.271267 -0.596847
+454 0.277013 -0.5911
+455 0.275867 -0.59742
+456 0.398853 -0.476153
+457 0.4 -0.470407
+458 0.405173 -0.465807
+459 0.410347 -0.46294
+460 0.416093 -0.459487
+461 0.421267 -0.45662
+462 0.42644 -0.45374
+463 0.43276 -0.450873
+464 0.439653 -0.447993
+465 0.446547 -0.446273
+466 0.45172 -0.44282
+467 0.456893 -0.43938
+468 0.46264 -0.435927
+469 0.469533 -0.433047
+470 0.475867 -0.43018
+471 0.481613 -0.431327
+472 0.477013 -0.435927
+473 0.47184 -0.438807
+474 0.466667 -0.441673
+475 0.467813 -0.44742
+476 0.462067 -0.450873
+477 0.457467 -0.454313
+478 0.452867 -0.457767
+479 0.44196 -0.459487
+480 0.435627 -0.462367
+481 0.430453 -0.465233
+482 0.425293 -0.468113
+483 0.43276 -0.472127
+484 0.425867 -0.47386
+485 0.42012 -0.47558
+486 0.414373 -0.477873
+487 0.408053 -0.480753
+488 0.402293 -0.47846
+489 0.51552 -0.345127
+490 0.51552 -0.338807
+491 0.516667 -0.333047
+492 0.52184 -0.33018
+493 0.531613 -0.329033
+494 0.53908 -0.330753
+495 0.54712 -0.332473
+496 0.553453 -0.330753
+497 0.557467 -0.335927
+498 0.562067 -0.340527
+499 0.5546 -0.341673
+500 0.54712 -0.340527
+501 0.540227 -0.340527
+502 0.53276 -0.33938
+503 0.52644 -0.3411
+504 0.522987 -0.346847
+505 0.522987 -0.352593
+506 0.51724 -0.349713
+507 0.74196 -0.49914
+508 0.735053 -0.49742
+509 0.731613 -0.49282
+510 0.7362 -0.48938
+511 0.74196 -0.485927
+512 0.747707 -0.48478
+513 0.7546 -0.484207
+514 0.759773 -0.487647
+515 0.766667 -0.488807
+516 0.764373 -0.49454
+517 0.758627 -0.496847
+518 0.748853 -0.49914
+518 0
+1 1 2
+2 2 3
+3 3 4
+4 4 5
+5 5 6
+6 6 7
+7 7 8
+8 8 9
+9 9 10
+10 10 11
+11 11 12
+12 12 13
+13 13 14
+14 14 15
+15 15 16
+16 16 17
+17 17 18
+18 18 19
+19 19 20
+20 20 21
+21 21 22
+22 22 23
+23 23 24
+24 24 25
+25 25 26
+26 26 27
+27 27 28
+28 28 29
+29 29 30
+30 30 31
+31 31 32
+32 32 33
+33 33 34
+34 34 35
+35 35 36
+36 36 37
+37 37 38
+38 38 39
+39 39 40
+40 40 41
+41 41 42
+42 42 43
+43 43 44
+44 44 45
+45 45 46
+46 46 47
+47 47 48
+48 48 49
+49 49 50
+50 50 51
+51 51 52
+52 52 53
+53 53 54
+54 54 55
+55 55 56
+56 56 57
+57 57 58
+58 58 59
+59 59 60
+60 60 61
+61 61 62
+62 62 63
+63 63 64
+64 64 65
+65 65 66
+66 66 67
+67 67 68
+68 68 69
+69 69 70
+70 70 71
+71 71 72
+72 72 73
+73 73 74
+74 74 75
+75 75 76
+76 76 77
+77 77 78
+78 78 79
+79 79 80
+80 80 81
+81 81 82
+82 82 83
+83 83 84
+84 84 85
+85 85 86
+86 86 87
+87 87 88
+88 88 89
+89 89 90
+90 90 91
+91 91 92
+92 92 93
+93 93 94
+94 94 95
+95 95 96
+96 96 97
+97 97 98
+98 98 99
+99 99 100
+100 100 101
+101 101 102
+102 102 103
+103 103 104
+104 104 105
+105 105 106
+106 106 107
+107 107 108
+108 108 109
+109 109 110
+110 110 111
+111 111 112
+112 112 113
+113 113 114
+114 114 115
+115 115 116
+116 116 117
+117 117 118
+118 118 119
+119 119 120
+120 120 121
+121 121 122
+122 122 123
+123 123 124
+124 124 125
+125 125 126
+126 126 127
+127 127 128
+128 128 129
+129 129 130
+130 130 131
+131 131 132
+132 132 133
+133 133 134
+134 134 135
+135 135 136
+136 136 137
+137 137 138
+138 138 139
+139 139 140
+140 140 141
+141 141 142
+142 142 143
+143 143 144
+144 144 145
+145 145 146
+146 146 147
+147 147 148
+148 148 149
+149 149 150
+150 150 151
+151 151 152
+152 152 153
+153 153 154
+154 154 155
+155 155 156
+156 156 157
+157 157 158
+158 158 159
+159 159 160
+160 160 161
+161 161 162
+162 162 163
+163 163 164
+164 164 165
+165 165 166
+166 166 167
+167 167 168
+168 168 169
+169 169 170
+170 170 171
+171 171 172
+172 172 173
+173 173 174
+174 174 175
+175 175 176
+176 176 177
+177 177 178
+178 178 179
+179 179 180
+180 180 181
+181 181 182
+182 182 183
+183 183 184
+184 184 185
+185 185 186
+186 186 187
+187 187 188
+188 188 189
+189 189 190
+190 190 191
+191 191 192
+192 192 193
+193 193 194
+194 194 195
+195 195 196
+196 196 197
+197 197 198
+198 198 199
+199 199 200
+200 200 201
+201 201 202
+202 202 203
+203 203 204
+204 204 205
+205 205 206
+206 206 207
+207 207 208
+208 208 209
+209 209 210
+210 210 211
+211 211 212
+212 212 213
+213 213 214
+214 214 215
+215 215 216
+216 216 217
+217 217 218
+218 218 219
+219 219 220
+220 220 221
+221 221 222
+222 222 223
+223 223 224
+224 224 225
+225 225 226
+226 226 227
+227 227 228
+228 228 229
+229 229 230
+230 230 231
+231 231 232
+232 232 233
+233 233 234
+234 234 235
+235 235 236
+236 236 237
+237 237 238
+238 238 239
+239 239 240
+240 240 241
+241 241 242
+242 242 243
+243 243 244
+244 244 245
+245 245 246
+246 246 247
+247 247 248
+248 248 249
+249 249 250
+250 250 251
+251 251 252
+252 252 253
+253 253 254
+254 254 255
+255 255 256
+256 256 257
+257 257 258
+258 258 259
+259 259 260
+260 260 261
+261 261 262
+262 262 263
+263 263 264
+264 264 265
+265 265 266
+266 266 267
+267 267 268
+268 268 269
+269 269 270
+270 270 271
+271 271 272
+272 272 273
+273 273 274
+274 274 275
+275 275 276
+276 276 277
+277 277 278
+278 278 279
+279 279 280
+280 280 281
+281 281 282
+282 282 283
+283 283 284
+284 284 285
+285 285 286
+286 286 287
+287 287 288
+288 288 289
+289 289 290
+290 290 291
+291 291 292
+292 292 293
+293 293 294
+294 294 295
+295 295 296
+296 296 297
+297 297 298
+298 298 299
+299 299 300
+300 300 301
+301 301 302
+302 302 303
+303 303 304
+304 304 305
+305 305 306
+306 306 307
+307 307 308
+308 308 309
+309 309 310
+310 310 311
+311 311 312
+312 312 313
+313 313 314
+314 314 315
+315 315 316
+316 316 317
+317 317 318
+318 318 319
+319 319 320
+320 320 321
+321 321 322
+322 322 323
+323 323 324
+324 324 325
+325 325 326
+326 326 327
+327 327 328
+328 328 329
+329 329 330
+330 330 331
+331 331 332
+332 332 333
+333 333 334
+334 334 335
+335 335 336
+336 336 337
+337 337 338
+338 338 339
+339 339 340
+340 340 341
+341 341 342
+342 342 343
+343 343 344
+344 344 345
+345 345 346
+346 346 347
+347 347 348
+348 348 349
+349 349 350
+350 350 351
+351 351 352
+352 352 353
+353 353 354
+354 354 355
+355 355 356
+356 356 357
+357 357 358
+358 358 359
+359 359 360
+360 360 361
+361 361 362
+362 362 363
+363 363 364
+364 364 365
+365 365 366
+366 366 367
+367 367 368
+368 368 369
+369 369 370
+370 370 371
+371 371 372
+372 372 373
+373 373 374
+374 374 375
+375 375 376
+376 376 377
+377 377 378
+378 378 379
+379 379 380
+380 380 381
+381 381 382
+382 382 383
+383 383 384
+384 384 385
+385 385 386
+386 386 387
+387 387 388
+388 388 389
+389 389 390
+390 390 391
+391 391 392
+392 392 393
+393 393 394
+394 394 395
+395 395 396
+396 396 397
+397 397 398
+398 398 399
+399 399 400
+400 400 401
+401 401 402
+402 402 403
+403 403 404
+404 404 405
+405 405 406
+406 406 407
+407 407 408
+408 408 409
+409 409 410
+410 410 411
+411 411 412
+412 412 413
+413 413 414
+414 414 415
+415 415 416
+416 416 417
+417 417 418
+418 418 419
+419 419 420
+420 420 421
+421 421 422
+422 422 423
+423 423 424
+424 424 425
+425 425 426
+426 426 427
+427 427 428
+428 428 429
+429 429 430
+430 430 431
+431 431 432
+432 432 433
+433 433 434
+434 434 435
+435 435 1
+436 436 437
+437 437 438
+438 438 439
+439 439 440
+440 440 441
+441 441 442
+442 442 443
+443 443 444
+444 444 445
+445 445 446
+446 446 436
+447 447 448
+448 448 449
+449 449 450
+450 450 451
+451 451 447
+452 452 453
+453 453 454
+454 454 455
+455 455 452
+456 456 457
+457 457 458
+458 458 459
+459 459 460
+460 460 461
+461 461 462
+462 462 463
+463 463 464
+464 464 465
+465 465 466
+466 466 467
+467 467 468
+468 468 469
+469 469 470
+470 470 471
+471 471 472
+472 472 473
+473 473 474
+474 474 475
+475 475 476
+476 476 477
+477 477 478
+478 478 479
+479 479 480
+480 480 481
+481 481 482
+482 482 483
+483 483 484
+484 484 485
+485 485 486
+486 486 487
+487 487 488
+488 488 456
+489 489 490
+490 490 491
+491 491 492
+492 492 493
+493 493 494
+494 494 495
+495 495 496
+496 496 497
+497 497 498
+498 498 499
+499 499 500
+500 500 501
+501 501 502
+502 502 503
+503 503 504
+504 504 505
+505 505 506
+506 506 489
+507 507 508
+508 508 509
+509 509 510
+510 510 511
+511 511 512
+512 512 513
+513 513 514
+514 514 515
+515 515 516
+516 516 517
+517 517 518
+518 518 507
+6
+1 0.55 -0.335
+2 0.75 -0.49
+3 0.45 -0.45
+4 0.275 -0.595
+5 0.258 -0.613
+6 0.24 -0.635