1
// $Id: 2D_InitMesh.cpp,v 1.22 2006/01/06 00:34:25 geuzaine Exp $
3
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
5
// This program is free software; you can redistribute it and/or modify
6
// it under the terms of the GNU General Public License as published by
7
// the Free Software Foundation; either version 2 of the License, or
8
// (at your option) any later version.
10
// This program is distributed in the hope that it will be useful,
11
// but WITHOUT ANY WARRANTY; without even the implied warranty of
12
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
// GNU General Public License for more details.
15
// You should have received a copy of the GNU General Public License
16
// along with this program; if not, write to the Free Software
17
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
20
// Please report all bugs and problems to <gmsh@geuz.org>.
24
Generation du maillage initial 2D
26
Dans l'algorithme de maillage, on doit creer un maillage initial a l'aide
27
de contours orientes aire a doite. Deux types de problemes peuvent se presenter :
29
1) probleme courant : un triangle se trouve a l'exterieur du contour, il faut
30
que sa variable de localisation soit declaree externe.
32
2) probleme rare : une arete du contour n'a pas ete creee
34
On procede comme suit :
35
trouve le triangle, on procede de la facon suivante :
36
- on trie les aretes de tous les triangles ;
37
- Pour chaque contour on prend chaque arete orientee ;
38
- Si cette arete n'existe pas (bsearch return false) alors on cherche
39
les 2 triangles reliant cette arete et on les modifie
40
- Si cette arete possede 2 triangles adjacents, on declare externe
41
le triangle dont le troisieme noeud se trouve a gauche de l'arete.
46
probleme 1) + = point du contour
55
2 3 a gauche de 1-2 -> tr externe
58
probleme 2) arete inexistante 1-2
81
extern PointRecord *gPointArray;
83
static Tree_T *ETree = NULL, *EDToSwap = NULL;
84
static int pointA, pointB, Counter, Stagnant, StopNow;
97
void makepermut(int numpoints)
100
permut = (int *)Malloc(numpoints * sizeof(int)); //Memory leak! this is never freed...
101
for(i = 0; i < numpoints; i++) {
102
ip = gPointArray[i].permu;
109
double Xma, Xmb, Yma, Ymb, det, t, q;
110
double XMa, XMb, YMa, YMb;
111
double xa1, xa2, ya1, ya2;
112
double xb1, xb2, yb1, yb2;
114
xa1 = gPointArray[pointA].where.h;
115
xa2 = gPointArray[pointB].where.h;
116
ya1 = gPointArray[pointA].where.v;
117
ya2 = gPointArray[pointB].where.v;
119
xb1 = gPointArray[e->from].where.h;
120
xb2 = gPointArray[e->to].where.h;
121
yb1 = gPointArray[e->from].where.v;
122
yb2 = gPointArray[e->to].where.v;
124
Xma = DMIN(xa1, xa2);
125
Xmb = DMIN(xb1, xb2);
127
XMa = DMAX(xa1, xa2);
128
XMb = DMAX(xb1, xb2);
130
if(XMa < Xmb || XMb < Xma)
133
Yma = DMIN(ya1, ya2);
134
Ymb = DMIN(yb1, yb2);
136
YMa = DMAX(ya1, ya2);
137
YMb = DMAX(yb1, yb2);
139
if(YMa < Ymb || YMb < Yma)
142
det = (xa2 - xa1) * (yb1 - yb2) - (xb1 - xb2) * (ya2 - ya1);
148
t = ((xb1 - xa1) * (yb1 - yb2) - (xb1 - xb2) * (yb1 - ya1)) / det;
149
q = ((xa2 - xa1) * (yb1 - ya1) - (ya2 - ya1) * (xb1 - xa1)) / det;
156
//printf("t=%g q=%g det=%g\n", t, q, det);
161
int compareED2(const void *a, const void *b)
166
if(q->NbOccur > w->NbOccur)
168
if(q->NbOccur < w->NbOccur)
170
if(q->from > w->from)
172
if(q->from < w->from)
181
int compareED(const void *a, const void *b)
186
if(q->from > w->from)
188
if(q->from < w->from)
197
void DoWeSwapED(void *data, void *dummy)
203
if(e->from == pointA || e->from == pointB ||
204
e->to == pointA || e->to == pointB) {
207
else if(Tree_PQuery(EDToSwap, e)) {
209
else if(crossED(e)) {
210
Tree_Add(EDToSwap, e);
214
void Print(void *data, void *dummy)
219
Msg(INFO, "%d %d", e->from, e->to);
222
void SwapED(void *data, void *dummy)
225
int from = 0, to = 0;
229
if(Stagnant && Counter <= StopNow)
234
if(!e->Liste[0] || !e->Liste[1]) {
235
Msg(GERROR, "Initial mesh is wrong. Try the new isotropic algorithm.");
239
if(e->from != e->Liste[0]->t.a && e->from != e->Liste[0]->t.b &&
240
e->from != e->Liste[0]->t.c)
242
if(e->from != e->Liste[1]->t.a && e->from != e->Liste[1]->t.b &&
243
e->from != e->Liste[1]->t.c)
245
if(e->to != e->Liste[0]->t.a && e->to != e->Liste[0]->t.b &&
246
e->to != e->Liste[0]->t.c)
248
if(e->to != e->Liste[1]->t.a && e->to != e->Liste[1]->t.b &&
249
e->to != e->Liste[1]->t.c)
253
if(e->from != e->Liste[0]->t.a && e->to != e->Liste[0]->t.a)
254
from = e->Liste[0]->t.a;
255
else if(e->from != e->Liste[0]->t.b && e->to != e->Liste[0]->t.b)
256
from = e->Liste[0]->t.b;
257
else if(e->from != e->Liste[0]->t.c && e->to != e->Liste[0]->t.c)
258
from = e->Liste[0]->t.c;
260
if(e->from != e->Liste[1]->t.a && e->to != e->Liste[1]->t.a)
261
to = e->Liste[1]->t.a;
262
else if(e->from != e->Liste[1]->t.b && e->to != e->Liste[1]->t.b)
263
to = e->Liste[1]->t.b;
264
else if(e->from != e->Liste[1]->t.c && e->to != e->Liste[1]->t.c)
265
to = e->Liste[1]->t.c;
274
e->Liste[0]->t.a = from;
275
e->Liste[0]->t.b = to;
276
e->Liste[0]->t.c = e->from;
278
e->Liste[1]->t.a = from;
279
e->Liste[1]->t.b = to;
280
e->Liste[1]->t.c = e->to;
282
e->from = IMAX(from, to);
283
e->to = IMIN(from, to);
286
void SuppressInETree(void *data, void *dummy)
288
if(!Tree_Suppress(ETree, data))
289
Msg(WARNING, "Cannot suppress in ETree");
292
void AddInETree(void *data, void *dummy)
294
Tree_Add(ETree, data);
297
int verifie_cas_scabreux(int pa, int pb, ContourRecord ** ListContours,
302
ContourRecord *contour;
305
Edge.from = IMAX(pa, pb);
306
Edge.to = IMIN(pa, pb);
307
for(k = 0; k < Nc; k++) {
308
contour = ListContours[k];
309
for(i = -1; i < (contour->numpoints - 1); i++) {
311
a = permut[contour->oriented_points[0].permu];
312
b = permut[contour->oriented_points[contour->numpoints - 1].permu];
315
a = permut[contour->oriented_points[i].permu];
316
b = permut[contour->oriented_points[i + 1].permu];
321
if(a != pa && b != pa && a != pb && b != pb)
329
void verify_edges(List_T * ListDelaunay, ContourRecord ** ListContour,
335
int ok, max, i, k, c, a, b;
336
ContourRecord *contour;
346
ETree = Tree_Create(sizeof(Edge), compareED);
347
for(i = 0; i < List_Nbr(ListDelaunay); i++) {
349
del_P = *(Delaunay **) List_Pointer(ListDelaunay, i);
351
Edge.from = IMAX(del_P->t.a, del_P->t.b);
352
Edge.to = IMIN(del_P->t.a, del_P->t.b);
355
if((pEdge = (ED *) Tree_PQuery(ETree, &Edge))) {
356
pEdge->Liste[1] = del_P;
359
Edge.Liste[0] = del_P;
360
Edge.Liste[1] = NULL;
361
Tree_Add(ETree, &Edge);
364
Edge.from = IMAX(del_P->t.a, del_P->t.c);
365
Edge.to = IMIN(del_P->t.a, del_P->t.c);
367
if((pEdge = (ED *) Tree_PQuery(ETree, &Edge))) {
368
pEdge->Liste[1] = del_P;
371
Edge.Liste[0] = del_P;
372
Edge.Liste[1] = NULL;
373
Tree_Add(ETree, &Edge);
376
Edge.from = IMAX(del_P->t.c, del_P->t.b);
377
Edge.to = IMIN(del_P->t.c, del_P->t.b);
379
if((pEdge = (ED *) Tree_PQuery(ETree, &Edge))) {
380
pEdge->Liste[1] = del_P;
383
Edge.Liste[0] = del_P;
384
Edge.Liste[1] = NULL;
385
Tree_Add(ETree, &Edge);
390
for(k = 0; k < NumContours; k++) {
391
contour = ListContour[k];
392
for(i = -1; i < contour->numpoints - 1; i++) {
395
a = permut[contour->oriented_points[0].permu];
396
b = permut[contour->oriented_points[contour->numpoints - 1].permu];
399
a = permut[contour->oriented_points[i].permu];
400
b = permut[contour->oriented_points[i + 1].permu];
403
Edge.from = IMAX(a, b);
404
Edge.to = IMIN(a, b);
406
if(!Tree_Search(ETree, &Edge))
415
Msg(INFO, "Swapping (%d missing edges)", ok);
419
Tree_Delete(EDToSwap);
420
EDToSwap = Tree_Create(sizeof(ED), compareED2);
421
for(k = 0; k < NumContours; k++) {
422
contour = ListContour[k];
423
for(i = -1; i < contour->numpoints - 1; i++) {
426
a = permut[contour->oriented_points[0].permu];
427
b = permut[contour->oriented_points[contour->numpoints - 1].permu];
430
a = permut[contour->oriented_points[i].permu];
431
b = permut[contour->oriented_points[i + 1].permu];
437
Tree_Action(ETree, DoWeSwapED);
441
Msg(INFO, "Elimination (%d swaps)", Tree_Nbr(EDToSwap));
443
Tree_Action(EDToSwap, SuppressInETree);
444
Tree_Action(EDToSwap, SwapED);
445
Tree_Action(EDToSwap, AddInETree);
446
} while(Tree_Nbr(EDToSwap));
448
Tree_Delete(EDToSwap);
455
int comparepermu(const void *i, const void *j)
457
return (gPointArray[*(PointNumero *) i].permu -
458
gPointArray[*(PointNumero *) j].permu);
461
void verify_inside(Delaunay * ListDelaunay, int NumDelaunay)
464
PointNumero pt[3], compare, a, b;
465
int i, inside, cont[3];
467
for(i = 0; i < NumDelaunay; i++) {
469
pt[0] = ListDelaunay[i].t.a;
470
pt[1] = ListDelaunay[i].t.b;
471
pt[2] = ListDelaunay[i].t.c;
472
qsort(pt, 3, sizeof(PointNumero), comparepermu);
474
cont[0] = gPointArray[pt[0]].numcontour;
475
cont[1] = gPointArray[pt[1]].numcontour;
476
cont[2] = gPointArray[pt[2]].numcontour;
479
if((cont[1] != cont[2]) || (cont[0] != cont[1]) || (cont[0] != cont[2]))
489
if((Is_left_of(a, b, compare) && (inside == 0))) {
490
ListDelaunay[i].t.position = EXTERN;
added
removed
Mesh/2D_InitMesh.cpp
