1
// $Id: 2D_Mesh_Triangle.cpp,v 1.15 2006/01/29 22:53:41 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>.
27
#if !defined(HAVE_TRIANGLE)
29
int Mesh_Triangle(Surface * s)
31
Msg(GERROR, "Triangle is not compiled in this version of Gmsh");
37
#define ANSI_DECLARATORS
49
void AddInMesh(Surface * sur, int nbbound, Vertex ** vertexbound,
50
struct triangulateio *out)
56
//Msg(INFO, "Add in database...");
58
vtable = (Vertex **) Malloc(out->numberofpoints * sizeof(Vertex *));
60
for(i = 0; i < nbbound; i++)
61
vtable[i] = vertexbound[i];
64
for(i = nbbound; i < out->numberofpoints; i++) {
65
vtable[i] = Create_Vertex(++(THEM->MaxPointNum),
66
out->pointlist[i * 2],
67
out->pointlist[i * 2 + 1], 0.0,
68
out->pointattributelist[i], 0.0);
69
Tree_Add(sur->Vertices, &vtable[i]);
73
Free(out->pointattributelist);
75
for(i = 0; i < out->numberoftriangles; i++) {
76
s = Create_Simplex(vtable[out->trianglelist[i * out->numberofcorners + 0]],
77
vtable[out->trianglelist[i * out->numberofcorners + 1]],
78
vtable[out->trianglelist[i * out->numberofcorners + 2]],
81
Tree_Add(sur->Simplexes, &s);
85
Free(out->trianglelist);
86
Free(out->triangleattributelist);
88
//Msg(INFO, "...done");
91
// This is horrible...
93
void FindPointInHole(List_T * verts, REAL * x, REAL * y)
96
double a[3], b[3], c[3];
97
List_Read(verts, 0, &v1);
98
List_Read(verts, 1, &v2);
99
a[0] = v2->Pos.X - v1->Pos.X;
100
a[1] = v2->Pos.Y - v1->Pos.Y;
107
*x = 0.5 * (v1->Pos.X + v2->Pos.X) + 1.e-12 * CTX.lc * c[0];
108
*y = 0.5 * (v1->Pos.Y + v2->Pos.Y) + 1.e-12 * CTX.lc * c[1];
111
int Mesh_Triangle(Surface * s)
114
int i, j, k, l, NbPts = 0, first;
118
struct triangulateio in, mid, out;
120
for(i = 0; i < List_Nbr(s->Contours); i++) {
121
List_Read(s->Contours, i, &list);
122
NbPts += List_Nbr(list);
125
in.numberofpoints = NbPts;
126
in.pointlist = (REAL *) Malloc(in.numberofpoints * 2 * sizeof(REAL));
127
vtable = (Vertex **) Malloc(in.numberofpoints * sizeof(Vertex *));
128
in.numberofpointattributes = 1;
129
in.pointattributelist = (REAL *) Malloc(in.numberofpoints *
130
in.numberofpointattributes *
132
in.pointmarkerlist = NULL;
133
in.numberofsegments = NbPts;
134
in.segmentlist = (int *)Malloc(in.numberofsegments * 2 * sizeof(int));
135
in.segmentmarkerlist = (int *)Malloc(in.numberofsegments * sizeof(int));
136
in.numberofregions = 0;
137
in.regionlist = NULL;
141
double lc_max = -1.0;
142
for(i = 0; i < List_Nbr(s->Contours); i++) {
143
List_Read(s->Contours, i, &list);
145
for(j = 0; j < List_Nbr(list); j++) {
146
List_Read(list, j, &v);
147
in.pointlist[k] = v->Pos.X;
148
in.pointlist[k + 1] = v->Pos.Y;
149
in.pointattributelist[l] = v->lc;
150
if(v->lc > lc_max) lc_max = v->lc;
152
in.segmentlist[k] = l;
153
in.segmentlist[k + 1] = (j == List_Nbr(list) - 1) ? (first) : (l + 1);
154
in.segmentmarkerlist[l] = i;
160
if(List_Nbr(s->Contours) > 1) {
161
in.numberofholes = List_Nbr(s->Contours) - 1;
162
in.holelist = (REAL *) Malloc(in.numberofholes * 2 * sizeof(REAL));
163
for(i = 1; i < List_Nbr(s->Contours); i++) {
164
List_Read(s->Contours, i, &list);
165
FindPointInHole(list, &in.holelist[(i - 1) * 2],
166
&in.holelist[(i - 1) * 2 + 1]);
170
in.numberofholes = 0;
174
mid.pointlist = NULL;
175
mid.pointattributelist = NULL;
176
mid.pointmarkerlist = NULL;
177
mid.trianglelist = NULL;
178
mid.triangleattributelist = NULL;
179
mid.neighborlist = NULL;
180
mid.segmentlist = NULL;
181
mid.segmentmarkerlist = NULL;
183
mid.edgemarkerlist = NULL;
185
// triangulate the points with minimum angle > 20 deg, with no
186
// boundary breaking, and with an area constraint related to the
187
// maximum char. length allowed (this last constraint is to avoid an
188
// extremely coarse initial grid, on which the interpolation of the
189
// final char. lengths would be awful).
191
double a = lc_max*lc_max / 1.2; // FIXME: bof
193
sprintf(opts, "pqzYa%f", a);
195
strcpy(opts, "pqzY");
196
if(CTX.verbosity < 4)
198
triangulate(opts, &in, &mid, NULL);
201
Free(in.pointattributelist);
202
Free(in.pointmarkerlist);
204
Free(in.segmentlist);
205
Free(in.segmentmarkerlist);
208
if(CTX.mesh.initial_only == 2) {
209
AddInMesh(s, NbPts, vtable, &mid);
213
mid.trianglearealist =
214
(REAL *) Malloc(mid.numberoftriangles * sizeof(REAL));
215
for(i = 0; i < mid.numberoftriangles; i++) {
216
if(THEM->BackgroundMeshType == ONFILE) {
217
double xx = 0.0, yy = 0.0;
218
for(j = 0; j < mid.numberofcorners; j++) {
219
k = mid.trianglelist[i * mid.numberofcorners + j];
220
xx += mid.pointlist[2 * k];
221
yy += mid.pointlist[2 * k + 1];
223
xx /= mid.numberofcorners;
224
yy /= mid.numberofcorners;
226
v = Create_Vertex(-1, xx, yy, 0.0, 0.0, 0.0);
227
Calcule_Z_Plan(&v, &dum);
228
Projette_Inverse(&v, &dum);
229
val = BGMXYZ(v->Pos.X, v->Pos.Y, v->Pos.Z);
230
val = val * val / 1.2; // FIXME: bof
235
for(j = 0; j < mid.numberofcorners; j++) {
236
k = mid.trianglelist[i * mid.numberofcorners + j];
237
val += mid.pointattributelist[k];
239
val /= mid.numberofcorners;
240
val = val * val / 1.5; // FIXME: bof
242
mid.trianglearealist[i] = val;
245
out.pointlist = NULL;
246
out.pointattributelist = NULL;
247
out.trianglelist = NULL;
248
out.triangleattributelist = NULL;
250
// refine the triangulation according to the triangle area
253
//strcpy(opts, "praqzBPY");
254
strcpy(opts, CTX.mesh.triangle_options);
255
if(CTX.verbosity < 4)
257
triangulate(opts, &mid, &out, NULL);
259
// free all allocated arrays + those allocated by Triangle
262
Free(mid.pointattributelist);
263
Free(mid.pointmarkerlist);
264
Free(mid.trianglelist);
265
Free(mid.triangleattributelist);
266
Free(mid.trianglearealist);
267
Free(mid.segmentlist);
268
Free(mid.segmentmarkerlist);
270
AddInMesh(s, NbPts, vtable, &out);
276
#endif // !HAVE_TRIANGLE
added
removed
Mesh/2D_Mesh_Triangle.cpp
