1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 /* SUMMARY : ... */
18 /* LICENSE : LGPLv3 */
19 /* ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE */
20 /* AUTHORS : Pascal Frey */
21 /* E-MAIL : pascal.frey@sorbonne-universite.fr */
22
23 #ifdef __cplusplus
24 extern "C" {
25 #endif
26
27 #include "medit.h"
28 #include "extern.h"
29 #include "sproto.h"
30
31 typedef struct saddle {
32 double r1, r2;
33 float pt[3];
34 float vp1[2], vp2[2];
35 int k;
36 } Saddle;
37
38 #define EPSD (float)1.e-14
39 #define EPS3 (float)1.e-03
40 #define MAXPS (int)100
41 #define HSIZ (float)0.001
42
43 static int idir[5] = {0, 1, 2, 0, 1};
44
closedBall(pMesh mesh,int depart,ubyte i)45 int closedBall(pMesh mesh, int depart, ubyte i) {
46 int adj, iadr;
47 ubyte voy;
48
49 voy = idir[i + 1];
50 iadr = 3 * (depart - 1) + 1;
51 adj = mesh->adja[iadr + voy];
52
53 /* search triangles in ball */
54 while (adj && adj != depart) {
55 voy = mesh->voy[iadr + voy];
56 iadr = 3 * (adj - 1) + 1;
57 voy = idir[voy + 2];
58 adj = mesh->adja[iadr + voy];
59 }
60
61 return (adj == depart);
62 }
63
listCritPoint(pScene sc,pMesh mesh)64 GLuint listCritPoint(pScene sc, pMesh mesh) {
65 pTriangle pt;
66 pPoint p0, p1, p2;
67 pSolution s0, s1, s2;
68 pMaterial pm;
69 Saddle sad[MAXPS];
70 GLuint dlist;
71 double aire, ux, uy, vx, vy, dd, cb0[3], cb1[3], cb2[3], vv[3][2], bc[3];
72 double rgb[3], a0, a1, delta, rr1, rr2, aa;
73 float p[3];
74 int *adj, iadr, i, i1, i2, k, m, ncp, ps, ifilt;
75 ubyte tag;
76 static double hsv[3] = {0.0f, 1.0f, 0.80f};
77
78 if (!mesh->nbb || mesh->nfield != mesh->dim) return (0);
79
80 if (mesh->nt && !hashTria(mesh)) return (0);
81
82 if (egal(sc->iso.val[0], sc->iso.val[MAXISO - 1])) return (0);
83
84 if (ddebug) printf("find critical points\n");
85
86 /* build list */
87 ncp = 0;
88 ps = 0;
89 dlist = glGenLists(1);
90 glNewList(dlist, GL_COMPILE);
91 if (glGetError( )) return (0);
92
93 glPointSize(4.0);
94
95 ifilt = 0;
96 hsv[0] = sc->iso.col[0];
97 hsvrgb(hsv, rgb);
98
99 for (m = 0; m < sc->par.nbmat; m++) {
100 pm = &sc->material[m];
101 k = pm->depmat[LTria];
102 if (!k || pm->flag) continue;
103
104 while (k != 0) {
105 pt = &mesh->tria[k];
106 if (!pt->v[0]) {
107 k = pt->nxt;
108 continue;
109 }
110
111 p0 = &mesh->point[pt->v[0]];
112 p1 = &mesh->point[pt->v[1]];
113 p2 = &mesh->point[pt->v[2]];
114 s0 = &mesh->sol[pt->v[0]];
115 s1 = &mesh->sol[pt->v[1]];
116 s2 = &mesh->sol[pt->v[2]];
117
118 ux = p1->c[0] - p0->c[0];
119 uy = p1->c[1] - p0->c[1];
120 vx = p2->c[0] - p0->c[0];
121 vy = p2->c[1] - p0->c[1];
122
123 aire = ux * vy - uy * vx;
124 if (aire == 0.0) {
125 k = pt->nxt;
126 continue;
127 } else if (aire < 0.0) {
128 p1 = &mesh->point[pt->v[2]];
129 p2 = &mesh->point[pt->v[1]];
130 s1 = &mesh->sol[pt->v[2]];
131 s2 = &mesh->sol[pt->v[1]];
132 aire = -aire;
133 }
134
135 /* coef des barycentriques */
136 aire = 1.0 / aire;
137 cb0[0] = p1->c[1] - p2->c[1];
138 cb0[1] = -(p1->c[0] - p2->c[0]);
139 cb0[2] = p1->c[0] * p2->c[1] - p1->c[1] * p2->c[0];
140
141 cb1[0] = p2->c[1] - p0->c[1];
142 cb1[1] = -(p2->c[0] - p0->c[0]);
143 cb1[2] = p2->c[0] * p0->c[1] - p2->c[1] * p0->c[0];
144
145 cb2[0] = p0->c[1] - p1->c[1];
146 cb2[1] = -(p0->c[0] - p1->c[0]);
147 cb2[2] = p0->c[0] * p1->c[1] - p0->c[1] * p1->c[0];
148
149 for (i = 0; i < 3; i++) {
150 vv[i][0] = aire * (cb0[i] * s0->m[0] + cb1[i] * s1->m[0] + cb2[i] * s2->m[0]);
151 vv[i][1] = aire * (cb0[i] * s0->m[1] + cb1[i] * s1->m[1] + cb2[i] * s2->m[1]);
152 }
153
154 dd = vv[0][0] * vv[1][1] - vv[0][1] * vv[1][0];
155 if (fabs(dd) < EPSD) {
156 k = pt->nxt;
157 continue;
158 }
159
160 dd = 1.0 / dd;
161 p[0] = dd * (vv[1][0] * vv[2][1] - vv[2][0] * vv[1][1]);
162 p[1] = dd * (vv[0][1] * vv[2][0] - vv[0][0] * vv[2][1]);
163 p[2] = 0.0;
164
165 if (p[0] < mesh->xmin - mesh->xtra || p[0] > mesh->xmax - mesh->xtra ||
166 p[1] < mesh->ymin - mesh->ytra || p[1] > mesh->ymax - mesh->ytra) {
167 k = pt->nxt;
168 continue;
169 } else if (!inTria(mesh, k, p, bc)) {
170 k = pt->nxt;
171 continue;
172 }
173
174 /* filtering boundary points */
175 tag = 0;
176
177 for (i = 0; i < 3; i++) {
178 tag |= (bc[i] < EPS3) << i;
179 }
180
181 if (tag) {
182 iadr = 3 * (k - 1) + 1;
183 adj = &mesh->adja[iadr];
184 ifilt++;
185
186 switch (tag) {
187 case 1:
188 if (!adj[0]) {
189 k = pt->nxt;
190 continue;
191 }
192
193 break;
194 case 2:
195 if (!adj[1]) {
196 k = pt->nxt;
197 continue;
198 }
199
200 break;
201 case 4:
202 if (!adj[2]) {
203 k = pt->nxt;
204 continue;
205 }
206
207 break;
208
209 case 3:
210 if (!closedBall(mesh, k, 2)) {
211 k = pt->nxt;
212 continue;
213 }
214
215 break;
216 case 5:
217 if (!closedBall(mesh, k, 1)) {
218 k = pt->nxt;
219 continue;
220 }
221
222 break;
223 case 6:
224 if (!closedBall(mesh, k, 0)) {
225 k = pt->nxt;
226 continue;
227 }
228
229 break;
230 }
231 }
232
233 /* eigenvalues of jacobian */
234 a1 = -(vv[0][0] + vv[1][1]);
235 a0 = vv[0][0] * vv[1][1] - vv[0][1] * vv[1][0];
236 delta = a1 * a1 - 4 * a0;
237 i1 = i2 = 0;
238 if (delta >= 0.0) {
239 delta = sqrt(delta);
240 rr1 = 0.5 * (-a1 + delta);
241 rr2 = 0.5 * (-a1 - delta);
242 } else {
243 delta = sqrt(fabs(delta));
244 rr1 = rr2 = -0.5 * a1;
245 i1 = i2 = 0.5 * delta;
246 }
247
248 /* classification */
249 if (i1 && i2) {
250 glColor3f(0.0, 1.0, 0.5);
251 if (rr1 == 0.0f && rr2 == 0.0f)
252 output2(p[0], p[1], "Cp"); /* center */
253 else if (rr1 > 0.0f && rr2 > 0.0f)
254 output2(p[0], p[1], "Rf"); /* repelling focus */
255 else if (rr1 < 0.0f && rr2 < 0.0f)
256 output2(p[0], p[1], "Af"); /* attracting focus */
257 } else if (!i1 && !i2) {
258 glColor3f(1.0, 0.5, 0.0);
259 if (rr1 > 0.0f && rr2 > 0.0f) {
260 output2(p[0], p[1], "Rn"); /* repelling node */
261 } else if (rr1 < 0.0f && rr2 < 0.0f) {
262 output2(p[0], p[1], "An"); /* attracting node */
263 } else if (rr1 * rr2 < 0.0f) {
264 output2(p[0], p[1], "Sp"); /* Saddle point */
265 if (ddebug) printf(" saddle point %f %f\n", p[0] + mesh->xtra, p[1] + mesh->ytra);
266
267 if (ps < MAXPS - 5) {
268 ++ps;
269 sad[ps].pt[0] = p[0];
270 sad[ps].pt[1] = p[1];
271 sad[ps].pt[2] = 0.0f;
272
273 /* eigenvalues */
274 sad[ps].r1 = rr1;
275 sad[ps].r2 = rr2;
276
277 /* eigenvectors */
278 aa = vv[0][0] * vv[0][0];
279 dd = vv[1][1] * vv[1][1];
280 delta = sqrt(aa - 2.0 * vv[0][0] * vv[1][1] + dd + 4.0 * vv[0][1] * vv[1][0]);
281 sad[ps].vp1[0] = -0.5 * (-vv[0][0] + vv[1][1] - delta);
282 sad[ps].vp1[1] = vv[0][1];
283
284 sad[ps].vp2[0] = -0.5 * (-vv[0][0] + vv[1][1] + delta);
285 sad[ps].vp2[1] = vv[0][1];
286
287 sad[ps].k = k;
288 }
289 }
290 }
291
292 /* point color */
293 glBegin(GL_POINTS);
294 glColor3dv(rgb);
295 glVertex2fv(p);
296 glEnd( );
297 pt->cpt--;
298 ++ncp;
299
300 k = pt->nxt;
301 }
302 }
303
304 glPointSize(1.0);
305 glEndList( );
306
307 if (ncp) fprintf(stdout, " %d critical points identified (%d filtered)\n", ncp, ifilt);
308
309 return (dlist);
310
311 if (ps) {
312 time_t t;
313
314 fprintf(stdout, " Building streamline(s)");
315 fflush(stdout);
316 t = clock( );
317
318 if (!sc->slist) {
319 sc->stream = createStream(sc, mesh);
320 if (!sc->stream) return (dlist);
321 }
322
323 for (k = 1; k <= ps; k++) {
324 if (ddebug) printf(" eigenv1 %f %f\n", sad[k].vp1[0], sad[k].vp1[1]);
325
326 listSaddleStream(sc, mesh, sad[k].k, sad[k].pt, sad[k].vp1, sad[k].r1);
327 sad[k].vp1[0] = -sad[k].vp1[0];
328 sad[k].vp1[1] = -sad[k].vp1[1];
329 listSaddleStream(sc, mesh, sad[k].k, sad[k].pt, sad[k].vp1, sad[k].r1);
330
331 if (ddebug) printf(" eigenv2 %f %f\n", sad[k].vp2[0], sad[k].vp2[1]);
332
333 listSaddleStream(sc, mesh, sad[k].k, sad[k].pt, sad[k].vp2, sad[k].r2);
334 sad[k].vp2[0] = -sad[k].vp2[0];
335 sad[k].vp2[1] = -sad[k].vp2[1];
336 listSaddleStream(sc, mesh, sad[k].k, sad[k].pt, sad[k].vp2, sad[k].r2);
337 }
338
339 sc->isotyp |= S_STREAML;
340 fprintf(stdout, ": %d lines", ps);
341 t = clock( ) - t;
342 fprintf(stdout, " %6.2f sec.\n", t / (float)CLOCKS_PER_SEC);
343 }
344
345 return (dlist);
346 }
347
348 #ifdef __cplusplus
349 }
350 #endif
351