diff --git a/Mesh/3D_Bricks.cpp b/Mesh/3D_Bricks.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fec85ea54023d090fa2556efbe492565ef6207c2
--- /dev/null
+++ b/Mesh/3D_Bricks.cpp
@@ -0,0 +1,97 @@
+// $Id: 3D_Bricks.cpp,v 1.18 2006-01-28 18:44:19 geuzaine Exp $
+//
+// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Gmsh.h"
+#include "Numeric.h"
+#include "Mesh.h"
+
+void AddSimplexInGrid(Mesh * m, Simplex * s, int boule_boite)
+{
+ double XminBox = 0., XmaxBox = 0., YminBox = 0., YmaxBox = 0., ZmaxBox = 0., ZminBox = 0.;
+ int Ix1, Ix2, Iy1, Iy2, Iz1, Iz2;
+ int i, j, k, index;
+ Brick Br, *pBrick;
+
+ if(!m->Grid.init) {
+ m->Grid.Bricks =
+ List_Create(m->Grid.Nx * m->Grid.Ny * m->Grid.Nz, 10, sizeof(Brick));
+ for(i = 0; i < m->Grid.Nx * m->Grid.Ny * m->Grid.Nz; i++) {
+ Br.pT = List_Create(2, 2, sizeof(Simplex *));
+ Br.N = i + 1;
+ List_Add(m->Grid.Bricks, &Br);
+ }
+ m->Grid.init = 1;
+ }
+
+ if(boule_boite == BOITE) {
+ XminBox = XmaxBox = s->V[0]->Pos.X;
+ YminBox = YmaxBox = s->V[0]->Pos.Y;
+ ZminBox = ZmaxBox = s->V[0]->Pos.Z;
+ for(i = 1; i < 4; i++) {
+ XminBox = DMIN(XminBox, s->V[i]->Pos.X);
+ XmaxBox = DMAX(XmaxBox, s->V[i]->Pos.X);
+ YminBox = DMIN(YminBox, s->V[i]->Pos.Y);
+ YmaxBox = DMAX(YmaxBox, s->V[i]->Pos.Y);
+ ZminBox = DMIN(ZminBox, s->V[i]->Pos.Z);
+ ZmaxBox = DMAX(ZmaxBox, s->V[i]->Pos.Z);
+ }
+ }
+ else if(boule_boite == BOULE) {
+ XminBox = s->Center.X - s->Radius;
+ XmaxBox = s->Center.X + s->Radius;
+ YminBox = s->Center.Y - s->Radius;
+ YmaxBox = s->Center.Y + s->Radius;
+ ZminBox = s->Center.Z - s->Radius;
+ ZmaxBox = s->Center.Z + s->Radius;
+ }
+
+
+ Ix1 = (int)((double)m->Grid.Nx * (XminBox - m->Grid.min.X) /
+ (m->Grid.max.X - m->Grid.min.X));
+ Ix2 = (int)((double)m->Grid.Nx * (XmaxBox - m->Grid.min.X) /
+ (m->Grid.max.X - m->Grid.min.X));
+ Iy1 = (int)((double)m->Grid.Ny * (YminBox - m->Grid.min.Y) /
+ (m->Grid.max.Y - m->Grid.min.Y));
+ Iy2 = (int)((double)m->Grid.Ny * (YmaxBox - m->Grid.min.Y) /
+ (m->Grid.max.Y - m->Grid.min.Y));
+ Iz1 = (int)((double)m->Grid.Nz * (ZminBox - m->Grid.min.Z) /
+ (m->Grid.max.Z - m->Grid.min.Z));
+ Iz2 = (int)((double)m->Grid.Nz * (ZmaxBox - m->Grid.min.Z) /
+ (m->Grid.max.Z - m->Grid.min.Z));
+
+ Ix1 = IMAX(Ix1, 0);
+ Ix2 = IMIN(Ix2, m->Grid.Nx - 1);
+ Iy1 = IMAX(Iy1, 0);
+ Iy2 = IMIN(Iy2, m->Grid.Ny - 1);
+ Iz1 = IMAX(Iz1, 0);
+ Iz2 = IMIN(Iz2, m->Grid.Nz - 1);
+
+ for(i = Ix1; i <= Ix2; i++) {
+ for(j = Iy1; j <= Iy2; j++) {
+ for(k = Iz1; k <= Iz2; k++) {
+ index = i + j * m->Grid.Nx + k * m->Grid.Nx * m->Grid.Ny;
+ pBrick = (Brick *) List_Pointer(m->Grid.Bricks, index);
+ List_Add(pBrick->pT, &s);
+ }
+ }
+ }
+
+}
diff --git a/Mesh/3D_Mesh_Old.cpp b/Mesh/3D_Mesh_Old.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3066460bce67bee795a8272374a55507891832f2
--- /dev/null
+++ b/Mesh/3D_Mesh_Old.cpp
@@ -0,0 +1,874 @@
+// $Id: 3D_Mesh_Old.cpp,v 1.16 2006-01-28 18:44:19 geuzaine Exp $
+//
+// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+/*
+ Isotropic Delaunay 3D
+
+ tant que l'arbre des tetraedres de qualites inacceptables
+ n'est pas vide {
+ prendre le plus mauvais tetraedre;
+ creer un nouveau point;
+ eliminer les tetraedres dont le cercle circonscrit contient le point;
+ reconstruire le volume convexe;
+ }
+
+*/
+
+
+#include "Gmsh.h"
+#include "Numeric.h"
+#include "Mesh.h"
+#include "3D_Mesh.h"
+#include "Create.h"
+#include "Context.h"
+
+extern Mesh *THEM, *LOCAL;
+extern Context_T CTX;
+extern int FACE_DIMENSION;
+
+static Tree_T *Tsd, *Sim_Sur_Le_Bord, *POINTS_TREE;
+static List_T *Simplexes_Destroyed, *Simplexes_New, *Suppress;
+static List_T *LLL, *POINTS;
+static Simplex *THES;
+static Vertex *THEV;
+static Tree_T *SimXFac;
+static double volume, LC3D;
+static int ZONEELIMINEE, Methode = 0;
+
+Simplex MyNewBoundary;
+int Alerte_Point_Scabreux;
+
+void DebugSimplexe(Simplex * s)
+{
+ int i;
+
+ fprintf(stderr, "Simplexe %p = %d %d %d %d \n",
+ s, s->V[0]->Num, s->V[1]->Num, s->V[2]->Num, s->V[3]->Num);
+
+ for(i = 0; i < 4; i++) {
+ if(s->S[i] != &MyNewBoundary)
+ printf(" face : %d %d %d -> Simplexe %p\n",
+ s->F[i].V[0]->Num, s->F[i].V[1]->Num, s->F[i].V[2]->Num,
+ s->S[i]);
+ else
+ printf(" face : %d %d %d -> Simplexe Boundary\n",
+ s->F[i].V[0]->Num, s->F[i].V[1]->Num, s->F[i].V[2]->Num);
+ }
+}
+
+void VSIM(void *a, void *b)
+{
+ Simplex *S;
+
+ S = *(Simplex **) a;
+ if(S->V[3])
+ volume += fabs(S->Volume_Simplexe());
+}
+
+void add_points(void *a, void *b)
+{
+ Tree_Insert(POINTS_TREE, a);
+}
+
+void add_points_2(void *a, void *b)
+{
+ List_Add(POINTS, a);
+}
+
+
+double Interpole_lcTetraedre(Simplex * s, Vertex * v)
+{
+ double mat[3][3], rhs[3], sol[3], det;
+
+ s->matsimpl(mat);
+ rhs[0] = v->Pos.X - s->V[0]->Pos.X;
+ rhs[1] = v->Pos.Y - s->V[0]->Pos.Y;
+ rhs[2] = v->Pos.Z - s->V[0]->Pos.Z;
+
+ sys3x3(mat, rhs, sol, &det);
+ if(det == 0.0 ||
+ (1. - sol[0] - sol[1] - sol[2]) > 1. ||
+ (1. - sol[0] - sol[1] - sol[2]) < 0. ||
+ sol[0] > 1. ||
+ sol[1] > 1. ||
+ sol[2] > 1. || sol[0] < 0. || sol[1] < 0. || sol[2] < 0.) {
+ return DMAX(s->V[0]->lc,
+ DMAX(s->V[1]->lc, DMAX(s->V[2]->lc, s->V[3]->lc)));
+ //sol[0] = sol[1] = sol[2] = 0.25;
+ }
+
+ return (s->V[0]->lc * (1. - sol[0] - sol[1] - sol[2]) +
+ sol[0] * s->V[1]->lc + sol[1] * s->V[2]->lc + sol[2] * s->V[3]->lc);
+}
+
+Vertex *NewVertex(Simplex * s)
+{
+ Vertex *v;
+
+ v =
+ Create_Vertex(++THEM->MaxPointNum, s->Center.X, s->Center.Y, s->Center.Z,
+ 1., 0.0);
+ v->lc = Interpole_lcTetraedre(s, v);
+
+ return (v);
+}
+
+int Pt_In_Circum(Simplex * s, Vertex * v)
+{
+ double d1, d2, eps;
+
+ /* Determine si un point est dans le cercle circonscrit a un simplexe */
+
+ d1 = s->Radius;
+ d2 = sqrt(DSQR(v->Pos.X - s->Center.X) +
+ DSQR(v->Pos.Y - s->Center.Y) + DSQR(v->Pos.Z - s->Center.Z));
+
+ eps = fabs(d1 - d2) / (d1 + d2);
+
+ if(eps < 1.e-12) {
+ return (0); // return 1 ? 0 ?
+ }
+
+ if(d2 < d1)
+ return (1);
+
+ return (0);
+}
+
+void Action_First_Simplexes(void *a, void *b)
+{
+ Simplex **q;
+
+ if(!THES) {
+ q = (Simplex **) a;
+ if(Pt_In_Circum(*q, THEV)) {
+ THES = *q;
+ }
+ }
+}
+
+void LiS(void *a, void *b)
+{
+ int j, N;
+ SxF SXF, *pSXF;
+ Simplex **pS, *S;
+
+ pS = (Simplex **) a;
+ S = *pS;
+ N = (S->V[3]) ? 4 : 3;
+
+ for(j = 0; j < N; j++) {
+ SXF.F = S->F[j];
+ if((pSXF = (SxF *) Tree_PQuery(SimXFac, &SXF))) {
+ /* Creation du lien */
+ S->S[j] = pSXF->S;
+ pSXF->S->S[pSXF->NumFaceSimpl] = S;
+ }
+ else {
+ SXF.S = S;
+ SXF.NumFaceSimpl = j;
+ Tree_Add(SimXFac, &SXF);
+ }
+ }
+}
+
+void RzS(void *a, void *b)
+{
+ int j, N;
+ Simplex **pS, *S;
+ pS = (Simplex **) a;
+ S = *pS;
+
+ N = (S->V[3]) ? 4 : 3;
+
+ for(j = 0; j < N; j++) {
+ if((S->S[j]) == NULL) {
+ S->S[j] = &MyNewBoundary;
+ }
+ }
+}
+
+/* Cree les liens entre les simplexes, c.a.d recherche les voisins */
+
+void Link_Simplexes(List_T * Sim, Tree_T * Tim)
+{
+ Simplex *S;
+ int i;
+
+ SimXFac = Tree_Create(sizeof(SxF), compareSxF);
+ if(Sim) {
+ for(i = 0; i < List_Nbr(Sim); i++) {
+ List_Read(Sim, i, &S);
+ LiS(&S, NULL);
+ }
+ for(i = 0; i < List_Nbr(Sim); i++) {
+ List_Read(Sim, i, &S);
+ RzS(&S, NULL);
+ }
+ }
+ else {
+ Tree_Action(Tim, LiS);
+ Tree_Action(Tim, RzS);
+ }
+ Tree_Delete(SimXFac);
+}
+
+void Box_6_Tetraedron(List_T * P, Mesh * m)
+{
+#define FACT 1.1
+#define LOIN 0.2
+
+ int i, j;
+ static int pts[8][3] = { {0, 0, 0},
+ {1, 0, 0},
+ {1, 1, 0},
+ {0, 1, 0},
+ {0, 0, 1},
+ {1, 0, 1},
+ {1, 1, 1},
+ {0, 1, 1}
+ };
+ static int tet[6][4] = { {1, 5, 2, 4},
+ {2, 5, 6, 4},
+ {4, 5, 6, 8},
+ {6, 4, 8, 7},
+ {6, 4, 7, 3},
+ {2, 3, 4, 6}
+ };
+ double Xm, Ym, Zm, XM, YM, ZM, Xc, Yc, Zc;
+ Simplex *S, *ps;
+ Vertex *V, *v, *pv;
+ List_T *smp;
+
+ smp = List_Create(8, 1, sizeof(Simplex *));
+
+ V = (Vertex *) Malloc(8 * sizeof(Vertex));
+
+ for(i = 0; i < List_Nbr(P); i++) {
+ List_Read(P, i, &v);
+ if(!i) {
+ Xm = XM = v->Pos.X;
+ Ym = YM = v->Pos.Y;
+ Zm = ZM = v->Pos.Z;
+ }
+ else {
+ Xm = DMIN(Xm, v->Pos.X);
+ XM = DMAX(XM, v->Pos.X);
+ Ym = DMIN(Ym, v->Pos.Y);
+ YM = DMAX(YM, v->Pos.Y);
+ Zm = DMIN(Zm, v->Pos.Z);
+ ZM = DMAX(ZM, v->Pos.Z);
+ }
+ }
+ if(Xm == XM)
+ XM = Xm + 1.;
+ if(Ym == YM)
+ YM = Ym + 1.;
+ if(Zm == ZM)
+ ZM = Zm + 1.;
+
+ Xc = XM - Xm;
+ Yc = YM - Ym;
+ Zc = ZM - Zm;
+
+ /* initialisation de la grille */
+
+ m->Grid.init = 0;
+ m->Grid.min.X = Xm - LOIN * FACT * Xc;
+ m->Grid.min.Y = Ym - LOIN * FACT * Yc;
+ m->Grid.min.Z = Zm - LOIN * FACT * Zc;
+ m->Grid.max.X = XM + LOIN * FACT * Xc;
+ m->Grid.max.Y = YM + LOIN * FACT * Yc;
+ m->Grid.max.Z = ZM + LOIN * FACT * Zc;
+
+ m->Grid.Nx = m->Grid.Ny = m->Grid.Nz = 20;
+
+ /* Longueur Caracteristique */
+
+ LC3D = sqrt(Xc * Xc + Yc * Yc + Zc * Zc);
+
+ /* Points de la boite de 1 a 8
+
+ Z 8____________7
+ | /| /|
+ | / | / |
+ | / | / |
+ 5|/___|________/6 |
+ | 4|________|___|3
+ | / | /
+ | / Y | /
+ | / | /
+ |/____________|/___ X
+ 1 2
+
+ */
+
+ for(i = 0; i < 8; i++) {
+ if(pts[i][0])
+ V[i].Pos.X = Xm - LOIN * Xc;
+ else
+ V[i].Pos.X = XM + LOIN * Xc;
+
+ if(pts[i][1])
+ V[i].Pos.Y = Ym - LOIN * Yc;
+ else
+ V[i].Pos.Y = YM + LOIN * Yc;
+
+ if(pts[i][2])
+ V[i].Pos.Z = Zm - LOIN * Zc;
+ else
+ V[i].Pos.Z = ZM + LOIN * Zc;
+
+ V[i].Num = -(++THEM->MaxPointNum);
+ pv = &V[i];
+ pv->lc = 1.0;
+ pv->Mov = NULL;
+ Tree_Replace(m->Vertices, &pv);
+ }
+
+ /* 6 Tetraedres forment le maillage de la boite */
+
+ for(i = 0; i < 6; i++) {
+ S = Create_Simplex(&V[tet[i][0] - 1], &V[tet[i][1] - 1],
+ &V[tet[i][2] - 1], &V[tet[i][3] - 1]);
+ List_Add(smp, &S);
+ }
+
+ Link_Simplexes(smp, NULL);
+ for(i = 0; i < List_Nbr(smp); i++) {
+ List_Read(smp, i, &ps);
+ for(j = 0; j < 4; j++)
+ if(ps->S[j] == NULL)
+ ps->S[j] = &MyNewBoundary;
+ Tree_Replace(m->Simplexes, &ps);
+ }
+
+}
+
+
+void Fill_Sim_Des(void *a, void *b)
+{
+ Simplex **S;
+ S = (Simplex **) a;
+ if(Pt_In_Circum(*S, THEV))
+ List_Add(Simplexes_Destroyed, a);
+}
+
+void TStoLS(void *a, void *b)
+{
+ List_Add(Simplexes_Destroyed, a);
+}
+
+void TAtoLA(void *a, void *b)
+{
+ List_Add(Simplexes_New, a);
+}
+
+void CrSi(void *a, void *b)
+{
+ SxF *S;
+ Simplex *s;
+ S = (SxF *) a;
+ if(S->NumFaceSimpl == 1) {
+ s = Create_Simplex(THEV, S->F.V[0], S->F.V[1], S->F.V[2]);
+ s->iEnt = ZONEELIMINEE;
+ THEM->Metric->setSimplexQuality(s);
+ List_Add(Simplexes_New, &s);
+ }
+ else if(S->NumFaceSimpl != 2) {
+ Msg(WARNING, "Huh! Panic in CrSi");
+ }
+}
+
+
+void NewSimplexes(Mesh * m, List_T * Sim, List_T * news)
+{
+ int i, j;
+ Tree_T *SimXFac;
+ Simplex *S;
+ SxF SXF, *pSXF;
+
+ SimXFac = Tree_Create(sizeof(SxF), compareSxF);
+
+ for(i = 0; i < List_Nbr(Sim); i++) {
+ List_Read(Sim, i, &S);
+ if(!i)
+ ZONEELIMINEE = S->iEnt;
+ else {
+ if(S->iEnt != ZONEELIMINEE) {
+ Msg(WARNING, "Huh! The elimination failed %d %d",
+ S->iEnt, ZONEELIMINEE);
+ }
+ }
+ for(j = 0; j < 4; j++) {
+ SXF.F = S->F[j];
+ if((pSXF = (SxF *) Tree_PQuery(SimXFac, &SXF))) {
+ (pSXF->NumFaceSimpl)++;
+ }
+ else {
+ SXF.NumFaceSimpl = 1;
+ Tree_Add(SimXFac, &SXF);
+ }
+ }
+ }
+
+ /* Les faces non communes sont obligatoirement a la frontiere ...
+ -> Nouveaux simplexes */
+
+ Tree_Action(SimXFac, CrSi);
+ Tree_Delete(SimXFac);
+}
+
+
+
+/* Methode recursive : Rempli Tsd les simplexes detruits
+ Invariant : Le simplexe est a eliminer
+ Le simplexe n'est pas encore considere */
+
+int recur_bowyer(Simplex * s)
+{
+ int i;
+
+ Tree_Insert(Tsd, &s);
+ for(i = 0; i < 4; i++) {
+ if(s->S[i] && s->S[i] != &MyNewBoundary && !Tree_Query(Tsd, &s->S[i])) {
+ if(Pt_In_Circum(s->S[i], THEV) && (s->iEnt == s->S[i]->iEnt)) {
+ recur_bowyer(s->S[i]);
+ }
+ else {
+ if(s->iEnt != s->S[i]->iEnt) {
+ //Msg(WARNING, "Point scabreux %d", s->S[i]->Num);
+ Alerte_Point_Scabreux = 1;
+ }
+ Tree_Insert(Sim_Sur_Le_Bord, &s->S[i]);
+ }
+ }
+ }
+ return 1;
+}
+
+bool Bowyer_Watson(Mesh * m, Vertex * v, Simplex * S, int force)
+{
+ int i;
+ Simplex *s;
+ double volumeold, volumenew;
+
+ THEV = v;
+
+ double x =
+ (S->V[0]->Pos.X + S->V[1]->Pos.X + S->V[2]->Pos.X + S->V[3]->Pos.X) / 4.;
+ double y =
+ (S->V[0]->Pos.Y + S->V[1]->Pos.Y + S->V[2]->Pos.Y + S->V[3]->Pos.Y) / 4.;
+ double z =
+ (S->V[0]->Pos.Z + S->V[1]->Pos.Z + S->V[2]->Pos.Z + S->V[3]->Pos.Z) / 4.;
+
+ if(force)
+ THEM->Metric->Identity();
+ else
+ THEM->Metric->setMetric(x, y, z);
+
+ Tsd = Tree_Create(sizeof(Simplex *), compareSimplex);
+ Sim_Sur_Le_Bord = Tree_Create(sizeof(Simplex *), compareSimplex);
+ List_Reset(Simplexes_Destroyed);
+ List_Reset(Simplexes_New);
+
+ if(Methode) {
+ Tree_Action(m->Simplexes, Fill_Sim_Des);
+ }
+ else {
+ recur_bowyer(S);
+ }
+
+ Tree_Action(Tsd, TStoLS);
+ NewSimplexes(m, Simplexes_Destroyed, Simplexes_New);
+
+ /* calcul des volumes des simplexes crees */
+
+ if(Alerte_Point_Scabreux || !CTX.mesh.speed_max) {
+ volume = 0.0;
+ for(i = 0; i < List_Nbr(Simplexes_Destroyed); i++) {
+ VSIM(List_Pointer(Simplexes_Destroyed, i), NULL);
+ }
+ volumeold = volume;
+ volume = 0.0;
+ for(i = 0; i < List_Nbr(Simplexes_New); i++) {
+ VSIM(List_Pointer(Simplexes_New, i), NULL);
+ }
+ volumenew = volume;
+ }
+ else {
+ volumeold = 1.0;
+ volumenew = 1.0;
+ }
+
+ /* critere du volume */
+
+ if((fabs(volumeold - volumenew) / (volumeold + volumenew)) > 1.e-8) {
+ if(Tree_Suppress(m->Simplexes, &S)) {
+ S->Quality = 0.0;
+ Tree_Add(m->Simplexes, &S);
+ }
+ if(force) {
+ List_Reset(Simplexes_Destroyed);
+ List_Reset(Simplexes_New);
+ Tree_Delete(Sim_Sur_Le_Bord);
+ Tree_Delete(Tsd);
+ //printf(" Aie Aie Aie volume changed %g -> %g\n",volumeold,volumenew);
+ return false;
+ }
+ }
+ else {
+ Tree_Add(m->Vertices, &THEV);
+ for(i = 0; i < List_Nbr(Simplexes_New); i++) {
+ Tree_Add(m->Simplexes, List_Pointer(Simplexes_New, i));
+ }
+
+ /* Suppression des simplexes elimines */
+
+ for(i = 0; i < List_Nbr(Simplexes_Destroyed); i++) {
+ List_Read(Simplexes_Destroyed, i, &s);
+ if(!Tree_Suppress(m->Simplexes, &s))
+ Msg(GERROR, "Impossible to delete simplex");
+ Free_Simplex(&s, 0);
+ }
+
+ /* Creation des liens entre nouveaux simplexes */
+
+ Tree_Action(Sim_Sur_Le_Bord, TAtoLA);
+ Link_Simplexes(Simplexes_New, m->Simplexes);
+ }
+
+ Tree_Delete(Sim_Sur_Le_Bord);
+ Tree_Delete(Tsd);
+ return true;
+}
+
+void Convex_Hull_Mesh(List_T * Points, Mesh * m)
+{
+ int i, j, N, n;
+ int Nbr_OK = 0, Nbr_NOTOK = 0;
+
+ N = List_Nbr(Points);
+ n = IMAX(N / 20, 1);
+
+ Box_6_Tetraedron(Points, m);
+ // List_Sort (Points, comparePosition);
+
+ for(i = 0; i < N; i++) {
+ THES = NULL;
+ List_Read(Points, i, &THEV);
+
+ if(Simplexes_New)
+ for(j = 0; j < List_Nbr(Simplexes_New); j++) {
+ Action_First_Simplexes(List_Pointer(Simplexes_New, j), NULL);
+ }
+
+ if(!THES) {
+ Tree_Action(m->Simplexes, Action_First_Simplexes);
+ Nbr_OK++;
+ }
+ else {
+ Nbr_NOTOK++;
+ }
+ if(i % n == n - 1) {
+ volume = 0.0;
+ Tree_Action(m->Simplexes, VSIM);
+ Msg(STATUS3, "Nod=%d/%d Elm=%d", i + 1, N, Tree_Nbr(m->Simplexes));
+ Msg(STATUS1, "Vol=%g", volume);
+ }
+ if(!THES) {
+ Msg(WARNING, "Vertex (%g,%g,%g) in no simplex", THEV->Pos.X,
+ THEV->Pos.Y, THEV->Pos.Z);
+ THEV->Pos.X +=
+ CTX.mesh.rand_factor * LC3D * (1. -
+ (double)rand() / (double)RAND_MAX);
+ THEV->Pos.Y +=
+ CTX.mesh.rand_factor * LC3D * (1. -
+ (double)rand() / (double)RAND_MAX);
+ THEV->Pos.Z +=
+ CTX.mesh.rand_factor * LC3D * (1. -
+ (double)rand() / (double)RAND_MAX);
+ Tree_Action(m->Simplexes, Action_First_Simplexes);
+ }
+ bool ca_marche = Bowyer_Watson(m, THEV, THES, 1);
+ int count = 0;
+ while(!ca_marche) {
+ count++;
+ double dx =
+ CTX.mesh.rand_factor * LC3D * (1. -
+ (double)rand() / (double)RAND_MAX);
+ double dy =
+ CTX.mesh.rand_factor * LC3D * (1. -
+ (double)rand() / (double)RAND_MAX);
+ double dz =
+ CTX.mesh.rand_factor * LC3D * (1. -
+ (double)rand() / (double)RAND_MAX);
+ THEV->Pos.X += dx;
+ THEV->Pos.Y += dy;
+ THEV->Pos.Z += dz;
+ THES = NULL;
+ Tree_Action(m->Simplexes, Action_First_Simplexes);
+ ca_marche = Bowyer_Watson(m, THEV, THES, 1);
+ THEV->Pos.X -= dx;
+ THEV->Pos.Y -= dy;
+ THEV->Pos.Z -= dz;
+ if(count > 5) {
+ N++;
+ List_Add(POINTS, &THEV);
+ Msg(WARNING, "Unable to add point %d (will try it later)", THEV->Num);
+ break;
+ }
+ }
+ }
+}
+
+void suppress_vertex(void *data, void *dum)
+{
+ Vertex **pv;
+
+ pv = (Vertex **) data;
+ if((*pv)->Num < 0)
+ List_Add(Suppress, pv);
+}
+
+void suppress_simplex(void *data, void *dum)
+{
+ Simplex **pv;
+
+ pv = (Simplex **) data;
+ if((*pv)->iEnt == 0)
+ List_Add(Suppress, pv);
+}
+
+void add_in_bgm(void *a, void *b)
+{
+ Simplex **s, *S;
+
+ s = (Simplex **) a;
+ S = *s;
+ List_Add(LLL, S);
+}
+
+void Bgm_With_Points(Mesh * bgm)
+{
+ int i;
+ Simplex *s;
+
+ bgm->BGM.bgm = List_Create(Tree_Nbr(bgm->Simplexes), 10, sizeof(Simplex));
+ LLL = bgm->BGM.bgm;
+ Tree_Action(bgm->Simplexes, add_in_bgm);
+ for(i = 0; i < List_Nbr(LLL); i++) {
+ s = (Simplex *) List_Pointer(LLL, i);
+ AddSimplexInGrid(bgm, s, BOITE);
+ }
+}
+
+void Create_BgMesh(int Type, double lc, Mesh * m)
+{
+ m->BGM.Typ = Type;
+ switch (Type) {
+ case CONSTANT:
+ m->BGM.lc = lc;
+ break;
+ case ONFILE:
+ break;
+ case WITHPOINTS:
+ m->BGM.bgm = NULL;
+ break;
+ }
+}
+
+void Maillage_Volume(void *data, void *dum)
+{
+ Volume *v, **pv;
+ Mesh M;
+ Surface S, *s;
+ Simplex *simp;
+ Vertex *newv;
+ int n, N;
+ double uvw[3];
+ int i;
+
+ FACE_DIMENSION = 2;
+
+ pv = (Volume **) data;
+ v = *pv;
+
+ if(v->Dirty) {
+ Msg(INFO, "Not meshing dirty Volume %d", v->Num);
+ return;
+ }
+
+ if(Extrude_Mesh(v)) {
+ }
+ else if(MeshTransfiniteVolume(v)) {
+ }
+ else if(v->Typ == 99999) {
+
+ Simplexes_New = List_Create(10, 10, sizeof(Simplex *));
+ Simplexes_Destroyed = List_Create(10, 10, sizeof(Simplex *));
+
+ LOCAL = &M;
+ Create_BgMesh(THEM->BGM.Typ, .2, LOCAL);
+ s = &S;
+
+ POINTS_TREE = Tree_Create(sizeof(Vertex *), comparePosition);
+ POINTS = List_Create(100, 100, sizeof(Vertex *));
+ LOCAL->Simplexes = v->Simplexes;
+ LOCAL->Vertices = v->Vertices;
+
+ for(i = 0; i < List_Nbr(v->Surfaces); i++) {
+ List_Read(v->Surfaces, i, &s);
+ Tree_Action(s->Vertices, add_points);
+ }
+ Tree_Action(POINTS_TREE, add_points_2);
+ Tree_Delete(POINTS_TREE);
+
+ N = List_Nbr(POINTS);
+ n = N / 30 + 1;
+
+ if(!N)
+ return;
+
+ /* Creation d'un maillage initial respectant la frontiere */
+
+ Msg(STATUS2, "Mesh 3D... (initial)");
+
+ Convex_Hull_Mesh(POINTS, LOCAL);
+
+ if(!Coherence(v, LOCAL))
+ Msg(GERROR,
+ "Surface recovery failed (send a bug report with the geo file to <gmsh@geuz.org>!)");
+
+ Link_Simplexes(NULL, LOCAL->Simplexes);
+
+ if(CTX.mesh.initial_only == 3) {
+ POINTS_TREE = THEM->Vertices;
+ Tree_Action(v->Vertices, add_points);
+ POINTS_TREE = THEM->Simplexes;
+ Tree_Action(v->Simplexes, add_points);
+ return;
+ }
+
+ /* Suppression des noeuds de num < 0 */
+
+ Suppress = List_Create(10, 10, sizeof(Vertex *));
+ Tree_Action(v->Vertices, suppress_vertex);
+ for(i = 0; i < List_Nbr(Suppress); i++) {
+ Tree_Suppress(v->Vertices, List_Pointer(Suppress, i));
+ }
+ List_Delete(Suppress);
+
+ /* Suppression des elements dont le num de vol == 0 (cad qui
+ n'appartiennent a auncun volume defini) */
+
+ Suppress = List_Create(10, 10, sizeof(Simplex *));
+ Tree_Action(v->Simplexes, suppress_simplex);
+ for(i = 0; i < List_Nbr(Suppress); i++) {
+ Tree_Suppress(v->Simplexes, List_Pointer(Suppress, i));
+ }
+
+ List_Delete(Suppress);
+
+ if(Tree_Nbr(LOCAL->Simplexes) == 0)
+ return;
+
+ /* Si il reste quelque chose a mailler en volume : */
+
+ Msg(STATUS2, "Mesh 3D... (final)");
+
+ v->Simplexes = LOCAL->Simplexes;
+
+ Bgm_With_Points(LOCAL);
+ POINTS_TREE = THEM->Simplexes;
+
+ Tree_Right(LOCAL->Simplexes, &simp);
+ i = 0;
+
+ while(simp->Quality > CONV_VALUE) {
+ newv = NewVertex(simp);
+ while(!simp->Pt_In_Simplexe(newv, uvw, 1.e-5) &&
+ (simp->S[0] == &MyNewBoundary ||
+ !simp->S[0]->Pt_In_Simplexe(newv, uvw, 1.e-5)) &&
+ (simp->S[1] == &MyNewBoundary ||
+ !simp->S[1]->Pt_In_Simplexe(newv, uvw, 1.e-5)) &&
+ (simp->S[2] == &MyNewBoundary ||
+ !simp->S[2]->Pt_In_Simplexe(newv, uvw, 1.e-5)) &&
+ (simp->S[3] == &MyNewBoundary ||
+ !simp->S[3]->Pt_In_Simplexe(newv, uvw, 1.e-5))) {
+ Tree_Suppress(LOCAL->Simplexes, &simp);
+ simp->Quality = 0.1;
+ Tree_Insert(LOCAL->Simplexes, &simp);
+ Tree_Right(LOCAL->Simplexes, &simp);
+ if(simp->Quality < CONV_VALUE)
+ break;
+ newv = NewVertex(simp);
+ }
+ if(simp->Quality < CONV_VALUE)
+ break;
+ i++;
+ if(i % n == n - 1) {
+ volume = 0.0;
+ Tree_Action(LOCAL->Simplexes, VSIM);
+ Msg(STATUS3, "Nod=%d Elm=%d",
+ Tree_Nbr(LOCAL->Vertices), Tree_Nbr(LOCAL->Simplexes));
+ Msg(STATUS1, "Vol(%g) Conv(%g->%g)", volume, simp->Quality,
+ CONV_VALUE);
+ double adv = 100. * (CONV_VALUE / simp->Quality);
+ }
+ Bowyer_Watson(LOCAL, newv, simp, 0);
+ Tree_Right(LOCAL->Simplexes, &simp);
+ }
+
+ POINTS_TREE = THEM->Vertices;
+ Tree_Action(v->Vertices, add_points);
+ POINTS_TREE = THEM->Simplexes;
+ Tree_Action(v->Simplexes, add_points);
+
+ if(CTX.mesh.quality) {
+ extern void SwapEdges3D(Mesh * M, Volume * v, double GammaPrescribed,
+ bool order);
+ Msg(STATUS3, "Swapping edges (1st pass)");
+ SwapEdges3D(THEM, v, CTX.mesh.quality, true);
+ Msg(STATUS3, "Swapping edges (2nd pass)");
+ SwapEdges3D(THEM, v, CTX.mesh.quality, false);
+ Msg(STATUS3, "Swapping edges (last pass)");
+ SwapEdges3D(THEM, v, CTX.mesh.quality, true);
+ }
+
+ if(CTX.mesh.nb_smoothing) {
+ /*
+ Msg(STATUS3, "Smoothing volume %d", v->Num);
+ tnxe = Tree_Create (sizeof (NXE), compareNXE);
+ create_NXE (v->Vertices, v->Simplexes, tnxe);
+ for (int i = 0; i < CTX.mesh.nb_smoothing; i++)
+ Tree_Action (tnxe, ActionLiss);
+ delete_NXE (tnxe);
+ Msg(STATUS3, "Swapping edges (last pass)");
+ SwapEdges3D (THEM, v, 0.5, true);
+ */
+ }
+
+ List_Delete(Simplexes_New);
+ List_Delete(Simplexes_Destroyed);
+ }
+
+}
diff --git a/Mesh/mGeomSearch.h b/Mesh/mGeomSearch.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d36f2c91794efb5c160befc64986328e9d9833b
--- /dev/null
+++ b/Mesh/mGeomSearch.h
@@ -0,0 +1,186 @@
+#ifndef _MGEOM_SEARCH_H
+#define _MGEOM_SEARCH_H
+
+#include <vector>
+#include <algorithm>
+
+#define MIN(x,y) ((x<y)?(x):(y))
+#define MAX(x,y) ((x>y)?(x):(y))
+#define TOL 1.e-06
+
+namespace AOMD {
+
+ template <class T>
+ class Brick
+ {
+ public:
+ std::vector<T*> Objects;
+ Brick(){}
+ T* operator [] (int i)
+ {
+ if(i < 0 || i >= Objects.size())throw i;
+ T *obj = Objects[i];
+ return obj;
+ }
+ int size(){return Objects.size();}
+ };
+
+ template <class T, class GetBox, class InteriorCheck>
+ class mGeomSearch
+ {
+ Brick<T> *theTable;
+ int Nx,Ny,Nz;
+ double Xmin,Xmax,Ymin,Ymax,Zmin,Zmax;
+ int getBrickId (double X, double Y, double Z)
+ {
+ int Ix = (int)((double)Nx * (X-Xmin) / (Xmax-Xmin));
+ int Iy = (int)((double)Ny * (Y-Ymin) / (Ymax-Ymin));
+ int Iz = (int)((double)Nz * (Z-Zmin) / (Zmax-Zmin));
+ Ix = MIN(Ix,Nx-1);
+ Iy = MIN(Iy,Ny-1);
+ Iz = MIN(Iz,Nz-1);
+ if(Ix<0)Ix=0;
+ if(Iy<0)Iy=0;
+ if(Iz<0)Iz=0;
+ int index = Ix + Iy * Nx + Iz * Nx * Ny;
+ return index;
+ }
+ Brick<T> * getBrick (int index)
+ {
+ if(index <0 || index >= Nx*Ny*Nz)throw index;
+ Brick<T> *b = &(theTable[index]);
+ return b;
+ }
+ GetBox getBox;
+ InteriorCheck interiorCheck;
+ public:
+ mGeomSearch (double x1,
+ double x2,
+ double y1,
+ double y2,
+ double z1,
+ double z2,
+ GetBox g,
+ InteriorCheck i,
+ int nx = 10,
+ int ny = 10,
+ int nz = 10) : Nx(nx), Ny(ny), Nz(nz),getBox(g),interiorCheck(i)
+ {
+ Xmin = x1-TOL; Xmax = x2+TOL;
+ Ymin = y1-TOL; Ymax = y2+TOL;
+ Zmin = z1-TOL; Zmax = z2+TOL;
+
+ theTable = new Brick<T> [Nx*Ny*Nz];
+ }
+
+ ~mGeomSearch ()
+ {
+ delete [] theTable;
+ };
+ Brick<T> *getBrick(double X, double Y, double Z)
+ {
+ return (getBrick(getBrickId(X,Y,Z)));
+ }
+ T * find (double X[3] , double U[3]);
+ bool AddObject ( T *);
+ bool RemoveObject ( T *);
+ };
+
+ template <class T,class GetBox, class InteriorCheck>
+ T * mGeomSearch<T,GetBox,InteriorCheck> :: find ( double X[3] ,
+ double U[3] )
+ {
+ Brick<T> *b = getBrick (X[0],X[1],X[2]);
+ if(!b) return 0;
+
+ for(int i=0;i<b->size();i++)
+ {
+ T * obj = (*b)[i];
+ if(interiorCheck(obj,X,U))return obj;
+ }
+ return 0;
+ }
+
+ template <class T,class GetBox, class InteriorCheck>
+ bool mGeomSearch<T,GetBox, InteriorCheck> :: AddObject ( T * obj )
+ {
+ int Ix1,Ix2,Iy1,Iy2,Iz1,Iz2;
+ int i,j,k,index;
+ Brick<T> *pBrick;
+
+ /*Template Objects must overload getBox function*/
+ double min[3];
+ double max[3];
+ getBox(obj,min,max);
+
+ Ix1 = (int)( (double)Nx * (min[0] - Xmin) /( Xmax - Xmin ));
+ Ix2 = (int)( (double)Nx * (max[0] - Xmin) /( Xmax - Xmin ));
+ Iy1 = (int)( (double)Ny * (min[1] - Ymin) /( Ymax - Ymin ));
+ Iy2 = (int)( (double)Ny * (max[1] - Ymin) /( Ymax - Ymin ));
+ Iz1 = (int)( (double)Nz * (min[2] - Zmin) /( Zmax - Zmin ));
+ Iz2 = (int)( (double)Nz * (max[2] - Zmin) /( Zmax - Zmin ));
+
+
+ Ix1 = MAX(Ix1,0);
+ Ix2 = MIN(Ix2,Nx-1);
+ Iy1 = MAX(Iy1,0);
+ Iy2 = MIN(Iy2,Ny-1);
+ Iz1 = MAX(Iz1,0);
+ Iz2 = MIN(Iz2,Nz-1);
+
+ for(i=Ix1;i<=Ix2;i++){
+ for(j=Iy1;j<=Iy2;j++){
+ for(k=Iz1;k<=Iz2;k++){
+ index = i + j * Nx + k * Nx * Ny;
+ pBrick = getBrick(index);
+ pBrick->Objects.push_back(obj);
+ }
+ }
+ }
+ return true;
+ }
+
+ template <class T,class GetBox, class InteriorCheck>
+ bool mGeomSearch<T,GetBox, InteriorCheck> :: RemoveObject ( T * obj )
+ {
+ int Ix1,Ix2,Iy1,Iy2,Iz1,Iz2;
+ int i,j,k,index;
+ Brick<T> *pBrick;
+
+ /*Template Objects must overload getBox function*/
+ double min[3];
+ double max[3];
+ getBox(obj,min,max);
+
+ Ix1 = (int)( (double)Nx * (min[0] - Xmin) /( Xmax - Xmin ));
+ Ix2 = (int)( (double)Nx * (max[0] - Xmin) /( Xmax - Xmin ));
+ Iy1 = (int)( (double)Ny * (min[1] - Ymin) /( Ymax - Ymin ));
+ Iy2 = (int)( (double)Ny * (max[1] - Ymin) /( Ymax - Ymin ));
+ Iz1 = (int)( (double)Nz * (min[2] - Zmin) /( Zmax - Zmin ));
+ Iz2 = (int)( (double)Nz * (max[2] - Zmin) /( Zmax - Zmin ));
+
+
+ Ix1 = MAX(Ix1,0);
+ Ix2 = MIN(Ix2,Nx-1);
+ Iy1 = MAX(Iy1,0);
+ Iy2 = MIN(Iy2,Ny-1);
+ Iz1 = MAX(Iz1,0);
+ Iz2 = MIN(Iz2,Nz-1);
+
+ for(i=Ix1;i<=Ix2;i++){
+ for(j=Iy1;j<=Iy2;j++){
+ for(k=Iz1;k<=Iz2;k++){
+ index = i + j * Nx + k * Nx * Ny;
+ pBrick = getBrick(index);
+ pBrick->Objects.erase ( std::remove (pBrick->Objects.begin(),pBrick->Objects.end(),obj) ,
+ pBrick->Objects.end () );
+ }
+ }
+ }
+ return true;
+ }
+
+} // end of namespace
+
+#undef TOL
+#endif
diff --git a/examples/cube.geo b/examples/cube.geo
deleted file mode 100644
index decdd4a696bb42dae6ae5bda42e91f3bb5b099e1..0000000000000000000000000000000000000000
--- a/examples/cube.geo
+++ /dev/null
@@ -1,83 +0,0 @@
-lcar1 = .0275;
-
-length = 0.25;
-height = 1.0;
-depth = 0.25;
-
-Point(newp) = {length/2,height/2,depth,lcar1}; /* Point 1 */
-Point(newp) = {length/2,height/2,0,lcar1}; /* Point 2 */
-Point(newp) = {-length/2,height/2,depth,lcar1}; /* Point 3 */
-Point(newp) = {-length/2,-height/2,depth,lcar1}; /* Point 4 */
-Point(newp) = {length/2,-height/2,depth,lcar1}; /* Point 5 */
-Point(newp) = {length/2,-height/2,0,lcar1}; /* Point 6 */
-Point(newp) = {-length/2,height/2,0,lcar1}; /* Point 7 */
-Point(newp) = {-length/2,-height/2,0,lcar1}; /* Point 8 */
-Line(1) = {3,1};
-Line(2) = {3,7};
-Line(3) = {7,2};
-Line(4) = {2,1};
-Line(5) = {1,5};
-Line(6) = {5,4};
-Line(7) = {4,8};
-Line(8) = {8,6};
-Line(9) = {6,5};
-Line(10) = {6,2};
-Line(11) = {3,4};
-Line(12) = {8,7};
-Line Loop(13) = {-6,-5,-1,11};
-Plane Surface(14) = {13};
-Line Loop(15) = {4,5,-9,10};
-Plane Surface(16) = {15};
-Line Loop(17) = {-3,-12,8,10};
-Plane Surface(18) = {17};
-Line Loop(19) = {7,12,-2,11};
-Plane Surface(20) = {19};
-Line Loop(21) = {-4,-3,-2,1};
-Plane Surface(22) = {21};
-Line Loop(23) = {8,9,6,7};
-Plane Surface(24) = {23};
-
-Surface Loop(25) = {14,24,-18,22,16,-20};
-Complex Volume(26) = {25};
-
-n = 21;
-
-/*
-sizex = 4;
-sizey = 6;
-sizez = 4;
-*/
-sizex = n*length;
-sizey = n*height;
-sizez = n*depth;
-
-sizex = 9;
-sizey = 33;
-sizez = 9;
-
-
-Transfinite Line {2,4,9,7} = sizez;
-Transfinite Line {11,5,10,12} = sizey;
-Transfinite Line {3,1,6,8} = sizex;
-
-Transfinite Surface {14} = {5,4,3,1};
-Transfinite Surface {16} = {6,2,1,5};
-Transfinite Surface {18} = {6,2,7,8};
-Transfinite Surface {20} = {3,7,8,4};
-Transfinite Surface {22} = {3,7,2,1};
-Transfinite Surface {24} = {4,8,6,5};
-
-Recombine Surface {14,16,18,20,22,24};
-
-Transfinite Volume {26} = {3,7,2,1,4,8,6,5};
-
-Physical Surface(200) = {14,16,18,20,24};
-
-Physical Volume(100) = {26};
-/*
-Mesh.use_cut_plane = 1;
-Mesh.cut_planea = 0;
-Mesh.cut_planeb = 0;
-Mesh.cut_planec = 1;
-Mesh.cut_planed = -0.125;
-*/
\ No newline at end of file