1 #include "medit.h"
2 #include "extern.h"
3 #include "sproto.h"
4
5 typedef struct saddle {
6 double r1,r2;
7 float pt[3];
8 float vp1[2],vp2[2];
9 int k;
10 } Saddle;
11
12
13 #define EPSD 1.e-14
14 #define EPS3 1.e-03
15 #define MAXPS 100
16 #define HSIZ 0.001
17
18 static int idir[5] = {0,1,2,0,1};
19
20
closedBall(pMesh mesh,int depart,ubyte i)21 int closedBall(pMesh mesh,int depart,ubyte i) {
22 pTriangle pt;
23 int adj,iadr;
24 ubyte voy;
25
26 pt = &mesh->tria[depart];
27 voy = idir[i+1];
28 iadr = 3*(depart-1)+1;
29 adj = mesh->adja[iadr+voy];
30
31 /* search triangles in ball */
32 while ( adj && adj != depart ) {
33 voy = mesh->voy[iadr+voy];
34 iadr = 3*(adj-1)+1;
35 voy = idir[voy+2];
36 adj = mesh->adja[iadr+voy];
37 }
38
39 return( adj == depart );
40 }
41
42
listCritPoint(pScene sc,pMesh mesh)43 GLuint listCritPoint(pScene sc,pMesh mesh) {
44 pTriangle pt;
45 pPoint p0,p1,p2;
46 pSolution s0,s1,s2;
47 pMaterial pm;
48 Saddle sad[MAXPS];
49 GLuint dlist;
50 double aire,ux,uy,vx,vy,dd,cb0[3],cb1[3],cb2[3],vv[3][2],bc[3];
51 double rgb[3],a0,a1,delta,rr1,rr2,aa,dmin;
52 float p[3];
53 int *adj,iadr,i,i1,i2,k,m,ncp,ps,ifilt;
54 ubyte typ,tag;
55 static double hsv[3] = {0.0f, 1.0f, 0.80f};
56 time_t t;
57
58 if ( !mesh->nbb || mesh->nfield != mesh->dim ) return(0);
59 if ( mesh->nt && !hashTria(mesh) ) return(0);
60 if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) ) return(0);
61 if ( ddebug ) printf("find critical points\n");
62
63 /* build list */
64 typ = 0;
65 ncp = 0;
66 ps = 0;
67 dlist = glGenLists(1);
68 glNewList(dlist,GL_COMPILE);
69 if ( glGetError() ) return(0);
70 glPointSize(4.0);
71
72 dmin = sc->dmax * EPS;
73 dmin *= dmin;
74 ifilt = 0;
75 hsv[0] = sc->iso.col[0];
76 hsvrgb(hsv,rgb);
77
78 for (m=0; m<sc->par.nbmat; m++) {
79 pm = &sc->material[m];
80 k = pm->depmat[LTria];
81 if ( !k || pm->flag ) continue;
82
83 while ( k != 0 ) {
84 pt = &mesh->tria[k];
85 if ( !pt->v[0] ) {
86 k = pt->nxt;
87 continue;
88 }
89
90 p0 = &mesh->point[pt->v[0]];
91 p1 = &mesh->point[pt->v[1]];
92 p2 = &mesh->point[pt->v[2]];
93 s0 = &mesh->sol[pt->v[0]];
94 s1 = &mesh->sol[pt->v[1]];
95 s2 = &mesh->sol[pt->v[2]];
96
97 ux = p1->c[0] - p0->c[0];
98 uy = p1->c[1] - p0->c[1];
99 vx = p2->c[0] - p0->c[0];
100 vy = p2->c[1] - p0->c[1];
101
102 aire = ux*vy - uy*vx;
103 if ( aire == 0.0 ) {
104 k = pt->nxt;
105 continue;
106 }
107 else if ( aire < 0.0 ) {
108 p1 = &mesh->point[pt->v[2]];
109 p2 = &mesh->point[pt->v[1]];
110 s1 = &mesh->sol[pt->v[2]];
111 s2 = &mesh->sol[pt->v[1]];
112 aire = -aire;
113 }
114
115 /* coef des barycentriques */
116 aire = 1.0 / aire;
117 cb0[0] = p1->c[1] - p2->c[1];
118 cb0[1] = -(p1->c[0] - p2->c[0]);
119 cb0[2] = p1->c[0]*p2->c[1] - p1->c[1]*p2->c[0];
120
121 cb1[0] = p2->c[1] - p0->c[1];
122 cb1[1] = -(p2->c[0] - p0->c[0]);
123 cb1[2] = p2->c[0]*p0->c[1] - p2->c[1]*p0->c[0];
124
125 cb2[0] = p0->c[1] - p1->c[1];
126 cb2[1] = -(p0->c[0] - p1->c[0]);
127 cb2[2] = p0->c[0]*p1->c[1] - p0->c[1]*p1->c[0];
128 for (i=0; i<3; i++) {
129 vv[i][0] = aire * (cb0[i]*s0->m[0] + cb1[i]*s1->m[0] + cb2[i]*s2->m[0]);
130 vv[i][1] = aire * (cb0[i]*s0->m[1] + cb1[i]*s1->m[1] + cb2[i]*s2->m[1]);
131 }
132
133 dd = vv[0][0]*vv[1][1] - vv[0][1]*vv[1][0];
134 if ( fabs(dd) < EPSD ) {
135 k = pt->nxt;
136 continue;
137 }
138
139 dd = 1.0 / dd;
140 p[0] = dd * (vv[1][0]*vv[2][1] - vv[2][0]*vv[1][1]);
141 p[1] = dd * (vv[0][1]*vv[2][0] - vv[0][0]*vv[2][1]);
142 p[2] = 0.0;
143
144 if ( p[0] < mesh->xmin-mesh->xtra || p[0] > mesh->xmax-mesh->xtra
145 || p[1] < mesh->ymin-mesh->ytra || p[1] > mesh->ymax-mesh->ytra ) {
146 k = pt->nxt;
147 continue;
148 }
149 else if ( !inTria(mesh,k,p,bc) ) {
150 k = pt->nxt;
151 continue;
152 }
153
154 /* filtering boundary points */
155 tag = 0;
156 for (i=0; i<3; i++) tag |= (bc[i] < EPS3) << i;
157 if ( tag ) {
158 iadr = 3*(k-1)+1;
159 adj = &mesh->adja[iadr];
160 ifilt ++;
161 switch (tag) {
162 case 1:
163 if ( !adj[0] ) {
164 k = pt->nxt;
165 continue;
166 }
167 break;
168 case 2:
169 if ( !adj[1] ) {
170 k = pt->nxt;
171 continue;
172 }
173 break;
174 case 4:
175 if ( !adj[2] ) {
176 k = pt->nxt;
177 continue;
178 }
179 break;
180
181 case 3:
182 if ( !closedBall(mesh,k,2) ) {
183 k = pt->nxt;
184 continue;
185 }
186 break;
187 case 5:
188 if ( !closedBall(mesh,k,1) ) {
189 k = pt->nxt;
190 continue;
191 }
192 break;
193 case 6:
194 if ( !closedBall(mesh,k,0) ) {
195 k = pt->nxt;
196 continue;
197 }
198 break;
199 }
200 }
201
202 /* eigenvalues of jacobian */
203 a1 = -(vv[0][0] + vv[1][1]);
204 a0 = vv[0][0] *vv[1][1] - vv[0][1]*vv[1][0];
205 delta = a1*a1 - 4*a0;
206 i1 = i2 = 0;
207 if ( delta >= 0.0 ) {
208 delta = sqrt(delta);
209 rr1 = 0.5 * (-a1+delta);
210 rr2 = 0.5 * (-a1-delta);
211 }
212 else {
213 delta = sqrt(fabs(delta));
214 rr1 = rr2 = -0.5 * a1;
215 i1 = i2 = 0.5 * delta;
216 }
217
218 /* classification */
219 if ( i1 && i2 ) {
220 glColor3f(0.0,1.0,0.5);
221 if ( rr1 == 0.0f && rr2 == 0.0f )
222 output2(p[0],p[1],"Cp"); /* center */
223 else if ( rr1 > 0.0f && rr2 > 0.0f )
224 output2(p[0],p[1],"Rf"); /* repelling focus */
225 else if ( rr1 < 0.0f && rr2 < 0.0f )
226 output2(p[0],p[1],"Af"); /* attracting focus */
227 }
228 else if ( !i1 && !i2 ) {
229 glColor3f(1.0,0.5,0.0);
230 if ( rr1 > 0.0f && rr2 > 0.0f )
231 output2(p[0],p[1],"Rn"); /* repelling node */
232 else if ( rr1 < 0.0f && rr2 < 0.0f )
233 output2(p[0],p[1],"An"); /* attracting node */
234 else if ( rr1*rr2 < 0.0f ) {
235 output2(p[0],p[1],"Sp"); /* Saddle point */
236 if ( ddebug )
237 printf(" saddle point %f %f\n",p[0]+mesh->xtra,p[1]+mesh->ytra);
238 if ( ps < MAXPS-5 ) {
239 ++ps;
240 sad[ps].pt[0] = p[0];
241 sad[ps].pt[1] = p[1];
242 sad[ps].pt[2] = 0.0f;
243
244 /* eigenvalues */
245 sad[ps].r1 = rr1;
246 sad[ps].r2 = rr2;
247
248 /* eigenvectors */
249 aa = vv[0][0]*vv[0][0];
250 dd = vv[1][1]*vv[1][1];
251 delta = sqrt(aa-2.0*vv[0][0]*vv[1][1]+dd+4.0*vv[0][1]*vv[1][0]);
252 sad[ps].vp1[0] = -0.5*(-vv[0][0]+vv[1][1]-delta);
253 sad[ps].vp1[1] = vv[0][1];
254
255 sad[ps].vp2[0] = -0.5*(-vv[0][0]+vv[1][1]+delta);
256 sad[ps].vp2[1] = vv[0][1];
257
258 sad[ps].k = k;
259 }
260 }
261 }
262
263 /* point color */
264 glBegin(GL_POINTS);
265 glColor3dv(rgb);
266 glVertex2fv(p);
267 glEnd();
268 pt->cpt--;
269 ++ncp;
270
271 k = pt->nxt;
272 }
273 }
274 glPointSize(1.0);
275 glEndList();
276
277 if ( ncp )
278 fprintf(stdout," %d critical points identified (%d filtered)\n",ncp,ifilt);
279
280 return(dlist);
281
282 if ( ps ) {
283 fprintf(stdout," Building streamline(s)");
284 fflush(stdout);
285 t = clock();
286
287 if ( !sc->slist ) {
288 sc->stream = createStream(sc,mesh);
289 if ( !sc->stream ) return(dlist);
290 }
291 for (k=1; k<=ps; k++) {
292 if ( ddebug ) printf(" eigenv1 %f %f\n",sad[k].vp1[0],sad[k].vp1[1]);
293
294 listSaddleStream(sc,mesh,sad[k].k,sad[k].pt,sad[k].vp1,sad[k].r1);
295 sad[k].vp1[0] = -sad[k].vp1[0];
296 sad[k].vp1[1] = -sad[k].vp1[1];
297 listSaddleStream(sc,mesh,sad[k].k,sad[k].pt,sad[k].vp1,sad[k].r1);
298
299 if ( ddebug ) printf(" eigenv2 %f %f\n",sad[k].vp2[0],sad[k].vp2[1]);
300 listSaddleStream(sc,mesh,sad[k].k,sad[k].pt,sad[k].vp2,sad[k].r2);
301 sad[k].vp2[0] = -sad[k].vp2[0];
302 sad[k].vp2[1] = -sad[k].vp2[1];
303 listSaddleStream(sc,mesh,sad[k].k,sad[k].pt,sad[k].vp2,sad[k].r2);
304 }
305 sc->isotyp |= S_STREAML;
306 fprintf(stdout,": %d lines",ps);
307 t = clock() - t;
308 fprintf(stdout," %6.2f sec.\n",t/(float)CLOCKS_PER_SEC);
309 }
310
311 return(dlist);
312 }
313