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