1 /*--------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* file: graphXO_2d.c */
7 /* */
8 /* */
9 /* description: simple graphical routines in 2d */
10 /* */
11 /*--------------------------------------------------------------------------*/
12 /* */
13 /* authors: Alfred Schmidt */
14 /* Zentrum fuer Technomathematik */
15 /* Fachbereich 3 Mathematik/Informatik */
16 /* Universitaet Bremen */
17 /* Bibliothekstr. 2 */
18 /* D-28359 Bremen, Germany */
19 /* */
20 /* Kunibert G. Siebert */
21 /* Institut fuer Mathematik */
22 /* Universitaet Augsburg */
23 /* Universitaetsstr. 14 */
24 /* D-86159 Augsburg, Germany */
25 /* */
26 /* http://www.mathematik.uni-freiburg.de/IAM/ALBERTA */
27 /* */
28 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
29 /* */
30 /*--------------------------------------------------------------------------*/
31
32 /****************************************************************************/
33 static const REAL_B vertex_lambda_2d[N_VERTICES_2D] = {
34 INIT_BARY_2D(1.0,0.0,0.0),
35 INIT_BARY_2D(0.0,1.0,0.0),
36 INIT_BARY_2D(0.0,0.0,1.0),
37 };
38
39 #if 0
40 static const REAL_B midedge_lambda[N_VERTICES_2D] = {
41 INIT_BARY_2D(0.0,0.5,0.5),
42 INIT_BARY_2D(0.5,0.0,0.5),
43 INIT_BARY_2D(0.5,0.5,0.0)
44 };
45 #endif
46
47 typedef struct min_max {
48 float xmin_2d[DIM_OF_WORLD], xmax_2d[DIM_OF_WORLD], diam_2d[DIM_OF_WORLD];
49 } MIN_MAX;
50
xminmax_fct_2d(const EL_INFO * elinfo,void * data)51 static void xminmax_fct_2d(const EL_INFO *elinfo, void *data)
52 {
53 MIN_MAX *ud = (MIN_MAX *)data;
54 int i, j;
55 PARAMETRIC *parametric = elinfo->mesh->parametric;
56
57 if (parametric)
58 {
59 REAL_D world[N_VERTICES_2D];
60 parametric->init_element(elinfo, parametric);
61 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D,
62 vertex_lambda_2d, world);
63
64 for (i = 0; i < N_VERTICES_2D; i++) {
65 for (j = 0; j < DIM_OF_WORLD; j++)
66 {
67 ud->xmin_2d[j] = MIN(ud->xmin_2d[j], world[i][j]);
68 ud->xmax_2d[j] = MAX(ud->xmax_2d[j], world[i][j]);
69 }
70 }
71 }
72 else {
73 DEBUG_TEST_FLAG(FILL_COORDS, elinfo);
74
75 for (i = 0; i < N_VERTICES_2D; i++)
76 for (j = 0; j < DIM_OF_WORLD; j++)
77 {
78 ud->xmin_2d[j] = MIN(ud->xmin_2d[j], elinfo->coord[i][j]);
79 ud->xmax_2d[j] = MAX(ud->xmax_2d[j], elinfo->coord[i][j]);
80 }
81 }
82 return;
83 }
84
85
graph_open_window_2d(const char * title,const char * geometry,REAL * world,MESH * mesh)86 static GRAPH_WINDOW graph_open_window_2d(const char *title,
87 const char *geometry,
88 REAL *world, MESH *mesh)
89 {
90 FUNCNAME("graph_open_window_2d");
91 MIN_MAX mm[1] = {{{0}}};
92 OGL_WINDOW *winO = NULL;
93 char Geometry[16];
94 int i;
95
96 if (world)
97 {
98 mm->xmin_2d[0] = world[0];
99 mm->xmax_2d[0] = world[1];
100 mm->xmin_2d[1] = world[2];
101 mm->xmax_2d[1] = world[3];
102 mm->diam_2d[0] = MAX(mm->xmax_2d[0] - mm->xmin_2d[0], 1.E-10);
103 mm->diam_2d[1] = MAX(mm->xmax_2d[1] - mm->xmin_2d[1], 1.E-10);
104 }
105 else if (mesh)
106 {
107 for (i=0; i<DIM_OF_WORLD; i++)
108 mm->xmax_2d[i] = -(mm->xmin_2d[i] = 1.0E10);
109 mesh_traverse(mesh, -1, CALL_LEAF_EL|FILL_COORDS, xminmax_fct_2d, mm);
110 for (i=0; i<DIM_OF_WORLD; i++)
111 {
112 mm->diam_2d[i] = MAX(mm->xmax_2d[i] - mm->xmin_2d[i], 1.E-10);
113 mm->xmin_2d[i] -= 0.1 * mm->diam_2d[i];
114 mm->xmax_2d[i] += 0.1 * mm->diam_2d[i];
115 mm->diam_2d[i] *= 1.2;
116 }
117 }
118 else
119 {
120 for (i=0; i<DIM_OF_WORLD; i++)
121 {
122 mm->xmin_2d[i] = 0.0;
123 mm->xmax_2d[i] = 1.0;
124 mm->diam_2d[i] = 1.0;
125 }
126 }
127
128 if (!title)
129 title = "ALBERTA graphics";
130
131 if (!geometry)
132 {
133 REAL SIZE = 400.0, xsize, ysize;
134
135 if (mm->diam_2d[0] >= mm->diam_2d[1])
136 {
137 xsize = SIZE;
138 ysize = SIZE * mm->diam_2d[1] / mm->diam_2d[0];
139 }
140 else
141 {
142 xsize = SIZE * mm->diam_2d[0] / mm->diam_2d[1];
143 ysize = SIZE;
144 }
145 snprintf(Geometry, 16, "%dx%d", (int)xsize, (int)ysize);
146 geometry = Geometry;
147 MSG("use geometry: %s\n", geometry);
148 }
149
150 winO = OGL_create_window(title, geometry);
151
152 TEST_EXIT(winO, "Could not create window!\n");
153
154 glMatrixMode(GL_PROJECTION);
155 glLoadIdentity();
156 glMatrixMode(GL_MODELVIEW);
157 glLoadIdentity();
158
159 for (i = 0; i < DIM_OF_WORLD; i++) {
160 winO->xmin[i] = mm->xmin_2d[i];
161 winO->xmax[i] = mm->xmax_2d[i];
162 }
163 glOrtho(winO->xmin[0], winO->xmax[0],
164 winO->xmin[1], winO->xmax[1],
165 -1.0, 1.0);
166
167 graph_clear_window((GRAPH_WINDOW) winO, rgb_white);
168
169 return((GRAPH_WINDOW) winO);
170 }
171
172 /****************************************************************************/
173
174 static FLAGS graph_mesh_flags_2d;
175 static const float *linecolor_2d;
176
177 #if HAVE_LIBGL
OGLgraph_mesh_fct_2d(const EL_INFO * elinfo,void * data)178 static void OGLgraph_mesh_fct_2d(const EL_INFO *elinfo, void *data)
179 {
180 EL *el = elinfo->el;
181 int i, j, i1,i2;
182 double xmid[DIM_OF_WORLD];
183 DOF **dof;
184 char cindex[20];
185 REAL_D world[N_VERTICES_2D];
186 const REAL_D *coord;
187 PARAMETRIC *parametric = elinfo->mesh->parametric;
188
189 if (parametric)
190 {
191 parametric->init_element(elinfo, parametric);
192 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D,
193 vertex_lambda_2d, world);
194 coord = (const REAL_D *)world;
195 }
196 else
197 {
198 DEBUG_TEST_FLAG(FILL_COORDS, elinfo);
199 coord = (const REAL_D *)elinfo->coord;
200 }
201
202 for (j=0; j<DIM_OF_WORLD;j++) {
203 xmid[j] = 0.0;
204 for (i=0;i<N_VERTICES_2D;i++) xmid[j] += coord[i][j];
205 xmid[j] /= N_VERTICES_2D;
206 }
207
208 /* show refinement/coarsening mark */
209 if (graph_mesh_flags_2d & GRAPH_MESH_ELEMENT_MARK)
210 {
211 if (el->mark > 0) glColor3fv(rgb_red);
212 else if (el->mark < 0) glColor3fv(rgb_blue);
213 else glColor3fv(rgb_white);
214
215 glBegin(GL_TRIANGLES);
216 glVertex2dv(coord[0]);
217 glVertex2dv(coord[1]);
218 glVertex2dv(coord[2]);
219 glEnd();
220 }
221
222
223 #if DEBUG
224 /* show element indices */
225 if (graph_mesh_flags_2d & GRAPH_MESH_ELEMENT_INDEX)
226 {
227 if (elinfo->level < 5) {
228 sprintf(cindex,"%d",el->index);
229 /* PLOT_text2d(xmid[0],xmid[1],cindex); */
230 }
231 }
232 #endif
233
234
235 if (graph_mesh_flags_2d & GRAPH_MESH_BOUNDARY)
236 {
237 /* show boundary edges */
238 DEBUG_TEST_FLAG(FILL_BOUND, elinfo);
239
240 for (i=0; i < N_WALLS_2D; i++) {
241 if (!IS_INTERIOR(wall_bound(elinfo, i))) {
242 if (linecolor_2d)
243 glColor3fv(linecolor_2d);
244 else if (IS_DIRICHLET(wall_bound(elinfo, i)))
245 glColor3fv(rgb_blue);
246 else
247 glColor3fv(rgb_red);
248
249 i1 = (i+1) % 3;
250 i2 = (i+2) % 3;
251
252 glBegin(GL_LINE_STRIP);
253 glVertex2dv(coord[i1]);
254 /* if (parametric) insert additional points... */
255 glVertex2dv(coord[i2]);
256 glEnd();
257 }
258 }
259 }
260 else
261 {
262 /* show all edges */
263 glColor3fv(linecolor_2d ? linecolor_2d : rgb_black);
264 glBegin(GL_LINE_LOOP);
265 for (i=0; i < N_VERTICES_2D; i++) {
266 glVertex2dv(coord[i]);
267 /* if (parametric) insert additional points... */
268 }
269 glEnd();
270 }
271
272 /* show dof[0] at vertices */
273 if (graph_mesh_flags_2d & GRAPH_MESH_VERTEX_DOF)
274 {
275 if ((dof = el->dof))
276 {
277 for (i=0; i<N_VERTICES_2D; i++)
278 {
279 if (dof[i]) sprintf(cindex,"%d",dof[i][0]);
280 }
281 }
282 }
283 return;
284 }
285 #endif
286
287
graph_mesh_2d(GRAPH_WINDOW win,MESH * mesh,const GRAPH_RGBCOLOR c,FLAGS flags)288 static void graph_mesh_2d(GRAPH_WINDOW win, MESH *mesh,
289 const GRAPH_RGBCOLOR c, FLAGS flags)
290 {
291 OGL_WINDOW *ogl_win = (OGL_WINDOW *) win;
292
293 if (!mesh) return;
294
295 linecolor_2d = c;
296 graph_mesh_flags_2d = flags;
297
298 OGL_set_std_window(ogl_win);
299
300 glLineWidth(1.0);
301
302 mesh_traverse(mesh, -1, CALL_LEAF_EL | FILL_COORDS | FILL_BOUND,
303 OGLgraph_mesh_fct_2d, NULL);
304 OGL_FLUSH(ogl_win);
305
306 return;
307 }
308
309 /****************************************************************************/
310
graph_close_window_2d(GRAPH_WINDOW win)311 static void graph_close_window_2d(GRAPH_WINDOW win)
312 {
313 OGL_destroy_window((OGL_WINDOW *)win);
314
315 return;
316 }
317
graph_clear_window_2d(GRAPH_WINDOW win,const GRAPH_RGBCOLOR c)318 static void graph_clear_window_2d(GRAPH_WINDOW win, const GRAPH_RGBCOLOR c)
319 {
320 OGL_clear_window((OGL_WINDOW *)win, c ? c : rgb_white);
321
322 return;
323 }
324
325 /****************************************************************************/
326 /****************************************************************************/
327
328 static REAL level_value = 0.0;
329 static int nrefine = 0;
330
331 static const DOF_REAL_VEC *drv = NULL;
332 static const DOF_REAL_D_VEC *drdv = NULL;
333
334 static const BAS_FCTS *bas_fcts = NULL;
335 static int n_bas_fcts = 0;
336 static const BAS_FCT *phi = NULL;
337 static const REAL *el_vec = NULL;
338 static const REAL_D *el_vec_d = NULL;
339 static PARAMETRIC *el_parametric = NULL;
340 static const EL_INFO *el_info = NULL;
341 static REAL val_min, val_max;
342
343 /****************************************************************************/
344 /* display of error estimators */
345 /****************************************************************************/
346
347 static REAL (*Get_el_est)(EL *el) = NULL;
348
el_est_fct(const EL_INFO * elinfo,const REAL * lambda)349 static REAL el_est_fct(const EL_INFO *elinfo, const REAL *lambda)
350 {
351 return(Get_el_est(elinfo->el));
352 }
353
graph_el_est_2d(GRAPH_WINDOW win,MESH * mesh,REAL (* get_el_est)(EL * el),REAL min,REAL max)354 static void graph_el_est_2d(GRAPH_WINDOW win, MESH *mesh,
355 REAL (*get_el_est)(EL *el),
356 REAL min, REAL max)
357 {
358 FUNCNAME("graph_el_est_2d");
359
360 TEST_EXIT(Get_el_est = get_el_est, "no get_el_est()\n");
361
362 graph_fvalues_2d(win, mesh, el_est_fct, 0, min, max, 0);
363 MSG("values in range [%.3le, %.3le]\n", val_min, val_max);
364 }
365
graph_level_recursive(int refine,const REAL * b[N_VERTICES_2D],REAL v[N_VERTICES_2D],const REAL * x[N_VERTICES_2D])366 static void graph_level_recursive(int refine, const REAL *b[N_VERTICES_2D],
367 REAL v[N_VERTICES_2D], const REAL *x[N_VERTICES_2D])
368 {
369 int i,j;
370
371 if (refine > 0) /* refine and call recursively */
372 {
373 const REAL *bnew[N_VERTICES_2D], *xnew[N_VERTICES_2D];
374 REAL_B bm[N_VERTICES_2D];
375 REAL vnew[N_VERTICES_2D], vm[N_VERTICES_2D], xm[N_VERTICES_2D][DIM_OF_WORLD];
376
377 for (j = 0; j < 3; j++)
378 {
379 bm[0][j] = 0.5*(b[1][j] + b[2][j]);
380 bm[1][j] = 0.5*(b[0][j] + b[2][j]);
381 bm[2][j] = 0.5*(b[0][j] + b[1][j]);
382 }
383
384 for (i = 0; i < 3; i++)
385 {
386 vm[i] = 0.0;
387 for (j=0; j<n_bas_fcts; j++)
388 vm[i] += el_vec[j] * PHI(bas_fcts, j, bm[i]); /* !!! phi[j](bm[i]); */
389 }
390
391 if (el_parametric)
392 {
393 el_parametric->coord_to_world(el_info, NULL, N_EDGES_2D,
394 (const REAL_B *)bm, xm);
395 }
396 else
397 {
398 for (j = 0; j < DIM_OF_WORLD; j++)
399 {
400 xm[0][j] = 0.5*(x[1][j] + x[2][j]);
401 xm[1][j] = 0.5*(x[0][j] + x[2][j]);
402 xm[2][j] = 0.5*(x[0][j] + x[1][j]);
403 }
404 }
405
406 bnew[0] = b[0]; bnew[1] = bm[2]; bnew[2] = bm[1];
407 xnew[0] = x[0]; xnew[1] = xm[2]; xnew[2] = xm[1];
408 vnew[0] = v[0]; vnew[1] = vm[2]; vnew[2] = vm[1];
409 graph_level_recursive(refine-1, bnew, vnew, xnew);
410
411 bnew[0] = bm[2]; bnew[1] = b[2]; bnew[2] = bm[1];
412 xnew[0] = xm[2]; xnew[1] = x[2]; xnew[2] = xm[1];
413 vnew[0] = vm[2]; vnew[1] = v[2]; vnew[2] = vm[1];
414 graph_level_recursive(refine-1, bnew, vnew, xnew);
415
416 bnew[0] = b[2]; bnew[1] = bm[2]; bnew[2] = bm[0];
417 xnew[0] = x[2]; xnew[1] = xm[2]; xnew[2] = xm[0];
418 vnew[0] = v[2]; vnew[1] = vm[2]; vnew[2] = vm[0];
419 graph_level_recursive(refine-1, bnew, vnew, xnew);
420
421 bnew[0] = bm[2]; bnew[1] = b[1]; bnew[2] = bm[0];
422 xnew[0] = xm[2]; xnew[1] = x[1]; xnew[2] = xm[0];
423 vnew[0] = vm[2]; vnew[1] = v[1]; vnew[2] = vm[0];
424 graph_level_recursive(refine-1, bnew, vnew, xnew);
425
426 return;
427 }
428
429 /* draw level line */
430
431 {
432 static REAL small = 1.0E-8;
433 int i1, i2, flag[3];
434 REAL s[3];
435 float xy[2];
436
437 for (i=0; i<3; i++) {
438 i1 = (i+1) % 3;
439 if (ABS(v[i1]-v[i]) >= small) {
440 s[i] = (level_value - v[i]) / (v[i1]-v[i]);
441 flag[i] = (s[i] <= 1.0 && s[i] >= 0.0);
442 }
443 else {
444 flag[i] = 0;
445 if (ABS(level_value - v[i]) <= small) {
446 glBegin(GL_LINE_STRIP);
447 xy[0] = x[i][0];
448 xy[1] = x[i][1];
449 glVertex2fv(xy);
450 xy[0] = x[i1][0];
451 xy[1] = x[i1][1];
452 glVertex2fv(xy);
453 glEnd();
454 }
455 }
456 }
457
458 for (i=0; i<3; i++) {
459 i1 = (i+1) % 3;
460 if (flag[i] && flag[i1]) {
461 i2 = (i+2) % 3;
462
463 glBegin(GL_LINE_STRIP);
464 for (j=0; j<2; j++)
465 xy[j] = x[i][j] + s[i] * (x[i1][j] - x[i][j]);
466 glVertex2fv(xy);
467 for (j=0; j<2; j++)
468 xy[j] = x[i1][j] + s[i1] * (x[i2][j] - x[i1][j]);
469 glVertex2fv(xy);
470 glEnd();
471 }
472 }
473 }
474 }
475
graph_level_fct(const EL_INFO * elinfo,void * data)476 static void graph_level_fct(const EL_INFO *elinfo, void *data)
477 {
478 /*FUNCNAME("graph_level_fct");*/
479 static REAL_B bnew[N_VERTICES_2D] = {
480 INIT_BARY_2D(1.0,0.0,0.0),
481 INIT_BARY_2D(0.0,1.0,0.0),
482 INIT_BARY_2D(0.0,0.0,1.0)
483 };
484
485 static const REAL *b[3]={bnew[0],bnew[1],bnew[2]};
486 const REAL *x[3];
487 REAL_D coord[N_VERTICES_2D];
488 REAL v[3];
489 EL *el = elinfo->el;
490 int i, j;
491 PARAMETRIC *parametric = elinfo->mesh->parametric;
492
493 el_vec = fill_el_real_vec(NULL, el, drv)->vec;
494
495 el_info = elinfo;
496 if (parametric)
497 {
498 el_parametric = parametric;
499 parametric->init_element(elinfo, parametric);
500 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D,
501 vertex_lambda_2d, coord);
502 for (i = 0; i < 3; i++)
503 x[i] = coord[i];
504 }
505 else
506 {
507 el_parametric = NULL;
508 for (i = 0; i < 3; i++)
509 x[i] = elinfo->coord[i];
510 }
511
512 for (i = 0; i < 3; i++)
513 {
514 v[i] = 0.0;
515 for (j=0; j<n_bas_fcts; j++)
516 v[i] += el_vec[j] * PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
517 }
518
519 graph_level_recursive(nrefine, b, v, x);
520 }
521
522 /****************************************************************************/
523
graph_level_2d(GRAPH_WINDOW win,const DOF_REAL_VEC * v,REAL level,const GRAPH_RGBCOLOR c,int refine)524 void graph_level_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v, REAL level,
525 const GRAPH_RGBCOLOR c, int refine)
526 {
527 FUNCNAME("graph_level_2d");
528 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
529
530 if (!v) return;
531
532 if (v->fe_space && v->fe_space->admin && v->fe_space->admin->mesh)
533 {
534 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
535 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
536 return;
537 }
538
539 OGL_set_std_window(ogl_win);
540
541 glLineWidth(1);
542 glColor3fv(c ? c : rgb_black);
543
544 bas_fcts = v->fe_space->bas_fcts;
545 n_bas_fcts = bas_fcts->n_bas_fcts;
546 phi = bas_fcts->phi;
547 drv = v;
548 level_value = level;
549 if (refine >= 0)
550 nrefine = refine;
551 else
552 nrefine = MAX(0, bas_fcts->degree-1);
553
554 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL | FILL_COORDS,
555 graph_level_fct, NULL);
556
557 OGL_FLUSH(ogl_win);
558 }
559 else
560 ERROR("no FE_SPACE OR DOF_ADMIN or MESH\n");
561 }
562
563 /****************************************************************************/
564
graph_level_d_recursive(int refine,const REAL * b[3],REAL v[3],const REAL * x[3])565 static void graph_level_d_recursive(int refine, const REAL *b[3], REAL v[3],
566 const REAL *x[3])
567 {
568 int i,j,k;
569
570 if (refine > 0) /* refine and call recursively */
571 {
572 const REAL *bnew[3], *xnew[3];
573 REAL_B bm[3];
574 REAL vnew[3], vm[3], xm[3][DIM_OF_WORLD], phi_b;
575 REAL_D vd;
576
577 for (j = 0; j < 3; j++)
578 {
579 bm[0][j] = 0.5*(b[1][j] + b[2][j]);
580 bm[1][j] = 0.5*(b[0][j] + b[2][j]);
581 bm[2][j] = 0.5*(b[0][j] + b[1][j]);
582 }
583
584 for (i = 0; i < 3; i++)
585 {
586 SET_DOW(0.0, vd);
587 for (j=0; j<n_bas_fcts; j++) {
588 phi_b = PHI(bas_fcts, j, bm[i]); /* !!! phi[j](bm[i]); */
589 for (k=0; k<DIM_OF_WORLD; ++k)
590 vd[k] += el_vec_d[j][k] * phi_b;
591 }
592 vm[i] = NORM_DOW(vd);
593 }
594
595 if (el_parametric)
596 {
597 el_parametric->coord_to_world(el_info, NULL, N_EDGES_2D,
598 (const REAL_B *)bm, xm);
599 }
600 else
601 {
602 for (j = 0; j < DIM_OF_WORLD; j++)
603 {
604 xm[0][j] = 0.5*(x[1][j] + x[2][j]);
605 xm[1][j] = 0.5*(x[0][j] + x[2][j]);
606 xm[2][j] = 0.5*(x[0][j] + x[1][j]);
607 }
608 }
609
610 bnew[0] = b[0]; bnew[1] = bm[2]; bnew[2] = bm[1];
611 xnew[0] = x[0]; xnew[1] = xm[2]; xnew[2] = xm[1];
612 vnew[0] = v[0]; vnew[1] = vm[2]; vnew[2] = vm[1];
613 graph_level_d_recursive(refine-1, bnew, vnew, xnew);
614
615 bnew[0] = bm[2]; bnew[1] = b[2]; bnew[2] = bm[1];
616 xnew[0] = xm[2]; xnew[1] = x[2]; xnew[2] = xm[1];
617 vnew[0] = vm[2]; vnew[1] = v[2]; vnew[2] = vm[1];
618 graph_level_d_recursive(refine-1, bnew, vnew, xnew);
619
620 bnew[0] = b[2]; bnew[1] = bm[2]; bnew[2] = bm[0];
621 xnew[0] = x[2]; xnew[1] = xm[2]; xnew[2] = xm[0];
622 vnew[0] = v[2]; vnew[1] = vm[2]; vnew[2] = vm[0];
623 graph_level_d_recursive(refine-1, bnew, vnew, xnew);
624
625 bnew[0] = bm[2]; bnew[1] = b[1]; bnew[2] = bm[0];
626 xnew[0] = xm[2]; xnew[1] = x[1]; xnew[2] = xm[0];
627 vnew[0] = vm[2]; vnew[1] = v[1]; vnew[2] = vm[0];
628 graph_level_d_recursive(refine-1, bnew, vnew, xnew);
629
630 return;
631 }
632
633 /* draw level line */
634
635 {
636 static REAL small = 1.0E-8;
637 int i1, i2, flag[3];
638 REAL s[3];
639 float xy[2];
640
641 for (i=0; i<3; i++) {
642 i1 = (i+1) % 3;
643 if (ABS(v[i1]-v[i]) >= small) {
644 s[i] = (level_value - v[i]) / (v[i1]-v[i]);
645 flag[i] = (s[i] <= 1.0 && s[i] >= 0.0);
646 }
647 else {
648 flag[i] = 0;
649 if (ABS(level_value - v[i]) <= small) {
650 glBegin(GL_LINE_STRIP);
651 xy[0] = x[i][0];
652 xy[1] = x[i][1];
653 glVertex2fv(xy);
654 xy[0] = x[i1][0];
655 xy[1] = x[i1][1];
656 glVertex2fv(xy);
657 glEnd();
658 }
659 }
660 }
661
662 for (i=0; i<3; i++) {
663 i1 = (i+1) % 3;
664 if (flag[i] && flag[i1]) {
665 i2 = (i+2) % 3;
666 glBegin(GL_LINE_STRIP);
667 for (j=0; j<2; j++)
668 xy[j] = x[i][j] + s[i] * (x[i1][j] - x[i][j]);
669 glVertex2fv(xy);
670 for (j=0; j<2; j++)
671 xy[j] = x[i1][j] + s[i1] * (x[i2][j] - x[i1][j]);
672 glVertex2fv(xy);
673 glEnd();
674 }
675 }
676 }
677 }
678
graph_level_d_fct(const EL_INFO * elinfo,void * data)679 static void graph_level_d_fct(const EL_INFO *elinfo, void *data)
680 {
681 FUNCNAME("graph_level_d_fct");
682 static REAL_B bnew[3]={INIT_BARY_2D(1.0,0.0,0.0),
683 INIT_BARY_2D(0.0,1.0,0.0),
684 INIT_BARY_2D(0.0,0.0,1.0)};
685 static const REAL *b[3]={bnew[0],bnew[1],bnew[2]};
686 const REAL *x[3];
687 REAL_D coord[N_VERTICES_2D], vd;
688 REAL v[3], phi_b;
689 EL *el = elinfo->el;
690 int i, j, k;
691 PARAMETRIC *parametric = elinfo->mesh->parametric;
692
693 if (bas_fcts->get_real_d_vec)
694 el_vec_d = fill_el_real_d_vec(NULL, el, drdv)->vec;
695 else
696 ERROR("no bas_fcts->get_real_d_vec()\n");
697
698 el_info = elinfo;
699 if (parametric)
700 {
701 el_parametric = parametric;
702 parametric->init_element(elinfo, parametric);
703 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D, vertex_lambda_2d,
704 coord);
705 for (i = 0; i < 3; i++)
706 x[i] = coord[i];
707 }
708 else {
709 el_parametric = NULL;
710 for (i = 0; i < 3; i++)
711 x[i] = elinfo->coord[i];
712 }
713
714 for (i = 0; i < 3; i++)
715 {
716 SET_DOW(0.0, vd);
717 for (j=0; j<n_bas_fcts; j++) {
718 phi_b = PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]);*/
719 for (k=0; k<DIM_OF_WORLD; ++k)
720 vd[k] += el_vec_d[j][k] * phi_b;
721 }
722 v[i] = NORM_DOW(vd);
723 }
724
725 graph_level_d_recursive(nrefine, b, v, x);
726 }
727
728 /****************************************************************************/
729
graph_level_d_2d(GRAPH_WINDOW win,const DOF_REAL_D_VEC * v,REAL level,const GRAPH_RGBCOLOR c,int refine)730 void graph_level_d_2d(GRAPH_WINDOW win, const DOF_REAL_D_VEC *v, REAL level,
731 const GRAPH_RGBCOLOR c, int refine)
732 {
733 FUNCNAME("graph_level_d_2d");
734 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
735
736 if (!v) return;
737
738 if (level < 0.0) return; /* no vector norms < 0 */
739
740 if (v->fe_space && v->fe_space->admin && v->fe_space->admin->mesh)
741 {
742 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
743 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
744 return;
745 }
746 OGL_set_std_window(ogl_win);
747 glLineWidth(1.0);
748 glColor3fv(c ? c : rgb_black);
749
750 bas_fcts = v->fe_space->bas_fcts;
751 n_bas_fcts = bas_fcts->n_bas_fcts;
752 phi = bas_fcts->phi;
753 drdv = v;
754 level_value = level;
755 if (refine >= 0)
756 nrefine = refine;
757 else
758 nrefine = MAX(0, bas_fcts->degree-1);
759
760 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL | FILL_COORDS,
761 graph_level_d_fct, NULL);
762
763 OGL_FLUSH(ogl_win);
764 }
765 else
766 ERROR("no FE_SPACE OR DOF_ADMIN or MESH\n");
767 }
768
769 /****************************************************************************/
val_minmax_fct(const EL_INFO * elinfo,void * data)770 static void val_minmax_fct(const EL_INFO *elinfo, void *data)
771 {
772 FUNCNAME("val_minmax_fct");
773 #define NTEST 4
774 static REAL_B b[NTEST] = {INIT_BARY_2D(1.0, 0.0, 0.0),
775 INIT_BARY_2D(0.0, 1.0, 0.0),
776 INIT_BARY_2D(0.0, 0.0, 1.0),
777 INIT_BARY_2D(0.333,0.333,0.334)};
778 int i, j, ntest;
779 REAL v;
780
781 if (bas_fcts->get_real_vec)
782 el_vec = fill_el_real_vec(NULL, elinfo->el, drv)->vec;
783 else
784 ERROR("no bas_fcts->get_real_vec()\n");
785
786 if (nrefine>0) ntest = NTEST;
787 else ntest = 3;
788
789 for (i = 0; i < ntest; i++)
790 {
791 v = 0.0;
792 for (j=0; j<n_bas_fcts; j++)
793 v += el_vec[j] * PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
794
795 val_min = MIN(val_min, v);
796 val_max = MAX(val_max, v);
797 }
798 #undef NTEST
799 }
800
val_d_minmax_fct(const EL_INFO * elinfo,void * data)801 static void val_d_minmax_fct(const EL_INFO *elinfo, void *data)
802 {
803 FUNCNAME("val_d_minmax_fct");
804 #define NTEST 4
805 static REAL_B b[NTEST] = {INIT_BARY_2D(1.0, 0.0, 0.0),
806 INIT_BARY_2D(0.0, 1.0, 0.0),
807 INIT_BARY_2D(0.0, 0.0, 1.0),
808 INIT_BARY_2D(0.333,0.333,0.334)};
809 int i, j, k, ntest;
810 REAL v, phi_b;
811 REAL_D vd;
812
813 if (bas_fcts->get_real_d_vec)
814 el_vec_d = fill_el_real_d_vec(NULL, elinfo->el, drdv)->vec;
815 else
816 ERROR("no bas_fcts->get_real_d_vec()\n");
817
818 if (nrefine>0) ntest = NTEST;
819 else ntest = 3;
820
821 for (i = 0; i < ntest; i++)
822 {
823 SET_DOW(0.0, vd);
824 for (j=0; j<n_bas_fcts; j++) {
825 phi_b = PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
826 for (k=0; k<DIM_OF_WORLD; ++k)
827 vd[k] += el_vec_d[j][k] * phi_b;
828 }
829 v = NORM_DOW(vd);
830
831 val_min = MIN(val_min, v);
832 val_max = MAX(val_max, v);
833 }
834 #undef NTEST
835 }
836
graph_levels_2d(GRAPH_WINDOW win,const DOF_REAL_VEC * v,int n,REAL const * levels,const GRAPH_RGBCOLOR * color,int refine)837 void graph_levels_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
838 int n, REAL const *levels, const GRAPH_RGBCOLOR *color,
839 int refine)
840 {
841 FUNCNAME("graph_levels_2d");
842 #define MAXLEV 100
843 REAL lev[MAXLEV];
844 GRAPH_RGBCOLOR col[MAXLEV];
845 const REAL *l;
846 const GRAPH_RGBCOLOR *c;
847 int nl, i;
848 float nlr;
849
850 if (!v) return;
851
852 nl = MIN(n, MAXLEV);
853 nlr = 1.0 / (float)(MAX(nl,1));
854
855 if (v->fe_space && v->fe_space->admin && v->fe_space->admin->mesh)
856 {
857 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
858 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
859 return;
860 }
861
862 if (levels)
863 {
864 l = levels;
865 }
866 else
867 {
868 drv = v;
869 bas_fcts = v->fe_space->bas_fcts;
870 n_bas_fcts = bas_fcts->n_bas_fcts;
871 phi = bas_fcts->phi;
872
873 if (refine >= 0)
874 nrefine = refine;
875 else
876 nrefine = MAX(0, bas_fcts->degree-1);
877
878 val_max = -(val_min = 1.0E20);
879 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL,
880 val_minmax_fct, NULL);
881 val_max = MAX(val_max, val_min + 1.0E-10);
882
883 for (i = 0; i < nl; i++)
884 lev[i] = val_min + ((i-0.5) * nlr) * (val_max-val_min);
885 l = lev;
886 }
887
888 if (color)
889 {
890 c = color;
891 }
892 else
893 {
894 c = (const GRAPH_RGBCOLOR *) col;
895 for (i = 0; i < nl; i++)
896 {
897 col[i][0] = i * nlr;
898 col[i][1] = 4.0 * (i * nlr) * (1.0 - i*nlr);
899 col[i][2] = 1.0 - i * nlr;
900 }
901 }
902
903 for (i = 0; i < nl; i++)
904 graph_level_2d(win, v, l[i], c[i], refine);
905 }
906 else
907 ERROR("no FE_SPACE or DOF_ADMIN or MESH\n");
908 }
909
910
graph_levels_d_2d(GRAPH_WINDOW win,const DOF_REAL_D_VEC * v,int n,REAL const * levels,const GRAPH_RGBCOLOR * color,int refine)911 void graph_levels_d_2d(GRAPH_WINDOW win, const DOF_REAL_D_VEC *v, int n,
912 REAL const *levels, const GRAPH_RGBCOLOR *color,
913 int refine)
914 {
915 FUNCNAME("graph_levels_d_2d");
916 #define MAXLEV 100
917 REAL lev[MAXLEV];
918 GRAPH_RGBCOLOR col[MAXLEV];
919 const GRAPH_RGBCOLOR *c;
920 const REAL *l;
921 int nl, i;
922 REAL nlr;
923
924 if (!v) return;
925
926 nl = MIN(n, MAXLEV);
927 nlr = 1.0 / (float)(MAX(nl,1));
928
929 if (v->fe_space && v->fe_space->admin && v->fe_space->admin->mesh)
930 {
931 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
932 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
933 return;
934 }
935
936 if (levels)
937 {
938 l = levels;
939 }
940 else
941 {
942 drdv = v;
943 bas_fcts = v->fe_space->bas_fcts;
944 n_bas_fcts = bas_fcts->n_bas_fcts;
945 phi = bas_fcts->phi;
946
947 if (refine >= 0)
948 nrefine = refine;
949 else
950 nrefine = MAX(0, bas_fcts->degree-1);
951
952 val_max = -(val_min = 1.0E20);
953 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL,
954 val_d_minmax_fct, NULL);
955 val_max = MAX(val_max, val_min + 1.0E-10);
956
957 for (i = 0; i < nl; i++)
958 lev[i] = val_min + ((i-0.5) * nlr) * (val_max-val_min);
959 l = lev;
960 }
961
962 if (color)
963 {
964 c = color;
965 }
966 else
967 {
968 c = (const GRAPH_RGBCOLOR *) col;
969 for (i = 0; i < nl; i++)
970 {
971 col[i][0] = i * nlr;
972 col[i][1] = 4.0 * (i * nlr) * (1.0 - i*nlr);
973 col[i][2] = 1.0 - i * nlr;
974 }
975 }
976
977 for (i = 0; i < nl; i++)
978 graph_level_d_2d(win, v, l[i], c[i], refine);
979 }
980 else
981 ERROR("no FE_SPACE or DOF_ADMIN or MESH\n");
982 }
983
984 /****************************************************************************/
985 static REAL val_fac = 1.0;
986
val_color(REAL v)987 static float *val_color(REAL v)
988 {
989 static GRAPH_RGBCOLOR c;
990 float vv = (float) ((v - val_min) * val_fac);
991
992 vv = MIN(1.0, vv);
993 vv = MAX(0.0, vv);
994
995 c[0] = vv;
996 c[1] = 4 * vv * (1.0 - vv);
997 c[2] = 1.0 - vv;
998
999 return(c);
1000 }
1001
1002
graph_value_recursive(int refine,const REAL * b[3],REAL v[3],const REAL * x[3])1003 static void graph_value_recursive(int refine, const REAL *b[3],
1004 REAL v[3], const REAL *x[3])
1005 {
1006 int i,j;
1007
1008 if (refine > 0) /* refine and call recursively */
1009 {
1010 const REAL *bnew[3], *xnew[3];
1011 REAL_B bm[3];
1012 REAL vnew[3], vm[3], xm[3][DIM_OF_WORLD];
1013
1014 for (j = 0; j < 3; j++)
1015 {
1016 bm[0][j] = 0.5*(b[1][j] + b[2][j]);
1017 bm[1][j] = 0.5*(b[0][j] + b[2][j]);
1018 bm[2][j] = 0.5*(b[0][j] + b[1][j]);
1019 }
1020
1021 for (i = 0; i < 3; i++)
1022 {
1023 for (vm[i] = j = 0; j < n_bas_fcts; j++)
1024 vm[i] += el_vec[j]*PHI(bas_fcts, j, bm[i]); /* !!! phi[j](bm[i]); */
1025 }
1026
1027 if (el_parametric)
1028 {
1029 el_parametric->coord_to_world(el_info, NULL, N_EDGES_2D,
1030 (const REAL_B *) bm, xm);
1031 }
1032 else
1033 {
1034 for (j = 0; j < DIM_OF_WORLD; j++)
1035 {
1036 xm[0][j] = 0.5*(x[1][j] + x[2][j]);
1037 xm[1][j] = 0.5*(x[0][j] + x[2][j]);
1038 xm[2][j] = 0.5*(x[0][j] + x[1][j]);
1039 }
1040 }
1041
1042 bnew[0] = b[0]; bnew[1] = bm[2]; bnew[2] = bm[1];
1043 xnew[0] = x[0]; xnew[1] = xm[2]; xnew[2] = xm[1];
1044 vnew[0] = v[0]; vnew[1] = vm[2]; vnew[2] = vm[1];
1045 graph_value_recursive(refine-1, bnew, vnew, xnew);
1046
1047 bnew[0] = bm[2]; bnew[1] = b[1]; bnew[2] = bm[0];
1048 xnew[0] = xm[2]; xnew[1] = x[1]; xnew[2] = xm[0];
1049 vnew[0] = vm[2]; vnew[1] = v[1]; vnew[2] = vm[0];
1050 graph_value_recursive(refine-1, bnew, vnew, xnew);
1051
1052 bnew[0] = b[2]; bnew[1] = bm[1]; bnew[2] = bm[0];
1053 xnew[0] = x[2]; xnew[1] = xm[1]; xnew[2] = xm[0];
1054 vnew[0] = v[2]; vnew[1] = vm[1]; vnew[2] = vm[0];
1055 graph_value_recursive(refine-1, bnew, vnew, xnew);
1056
1057 bnew[0] = bm[0]; bnew[1] = bm[1]; bnew[2] = bm[2];
1058 xnew[0] = xm[0]; xnew[1] = xm[1]; xnew[2] = xm[2];
1059 vnew[0] = vm[0]; vnew[1] = vm[1]; vnew[2] = vm[2];
1060 graph_value_recursive(refine-1, bnew, vnew, xnew);
1061
1062 return;
1063 }
1064
1065 /* draw triangle */
1066
1067 glBegin(GL_TRIANGLES);
1068 for (i = 0; i < N_VERTICES_2D; i++) {
1069 glColor3fv(val_color(v[i]));
1070 glVertex2dv(x[i]);
1071 }
1072
1073 glEnd();
1074 }
1075
1076
graph_value_fct(const EL_INFO * elinfo,void * data)1077 static void graph_value_fct(const EL_INFO *elinfo, void *data)
1078 {
1079 FUNCNAME("graph_value_fct");
1080 static REAL_B bnew[3]={INIT_BARY_2D(1.0,0.0,0.0),
1081 INIT_BARY_2D(0.0,1.0,0.0),
1082 INIT_BARY_2D(0.0,0.0,1.0)};
1083 static const REAL *b[3]={bnew[0],bnew[1],bnew[2]};
1084 static const REAL *x[3];
1085 static REAL_D coord[N_VERTICES_2D];
1086 static REAL v[3];
1087 int i, j;
1088 PARAMETRIC *parametric = elinfo->mesh->parametric;
1089
1090 if (bas_fcts->get_real_vec)
1091 el_vec = fill_el_real_vec(NULL, elinfo->el, drv)->vec;
1092 else
1093 ERROR("no bas_fcts->get_real_vec()\n");
1094
1095 el_info = elinfo;
1096 if (parametric)
1097 {
1098 el_parametric = parametric;
1099 parametric->init_element(elinfo, parametric);
1100 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D, vertex_lambda_2d,
1101 coord);
1102 for (i = 0; i < 3; i++)
1103 x[i] = coord[i];
1104 }
1105 else {
1106 el_parametric = NULL;
1107 for (i = 0; i < 3; i++)
1108 x[i] = elinfo->coord[i];
1109 }
1110
1111 for (i = 0; i < 3; i++)
1112 {
1113 v[i] = 0.0;
1114 for (j=0; j<n_bas_fcts; j++)
1115 v[i] += el_vec[j] * PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
1116 }
1117
1118 graph_value_recursive(nrefine, b, v, x);
1119 }
1120
1121
graph_drv_2d(GRAPH_WINDOW win,const DOF_REAL_VEC * v,REAL min,REAL max,int refine)1122 static void graph_drv_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
1123 REAL min, REAL max, int refine)
1124 {
1125 FUNCNAME("graph_drv_2d");
1126 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
1127
1128 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
1129 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
1130 return;
1131 }
1132
1133 drv = v;
1134 bas_fcts = v->fe_space->bas_fcts;
1135 n_bas_fcts = bas_fcts->n_bas_fcts;
1136 phi = bas_fcts->phi;
1137
1138 if (refine >= 0)
1139 nrefine = refine;
1140 else
1141 nrefine = MAX(0, bas_fcts->degree-1);
1142
1143 if (min < max) {
1144 val_min = min;
1145 val_max = max;
1146 }
1147 else {
1148 val_max = -(val_min = 1.0E20);
1149 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL,
1150 val_minmax_fct, NULL);
1151 val_max = MAX(val_max, val_min + 1.0E-10);
1152 }
1153 val_fac = 1.0 / (val_max - val_min);
1154
1155
1156 OGL_set_std_window(ogl_win);
1157
1158 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL|FILL_COORDS,
1159 graph_value_fct, NULL);
1160
1161 OGL_FLUSH(ogl_win);
1162 }
1163
1164 /****************************************************************************/
1165
graph_value_d_recursive(int refine,const REAL * b[3],REAL v[3],const REAL * x[3])1166 static void graph_value_d_recursive(int refine, const REAL *b[3],
1167 REAL v[3], const REAL *x[3])
1168 {
1169 int i, j, k;
1170 REAL phi_b;
1171
1172 if (refine > 0) /* refine and call recursively */
1173 {
1174 const REAL *bnew[3], *xnew[3];
1175 REAL_B bm[3];
1176 REAL vnew[3], vm[3], xm[3][DIM_OF_WORLD];
1177 REAL_D vd;
1178
1179 for (j = 0; j < 3; j++)
1180 {
1181 bm[0][j] = 0.5*(b[1][j] + b[2][j]);
1182 bm[1][j] = 0.5*(b[0][j] + b[2][j]);
1183 bm[2][j] = 0.5*(b[0][j] + b[1][j]);
1184 }
1185
1186 for (i = 0; i < 3; i++)
1187 {
1188 SET_DOW(0.0, vd);
1189 for (j = 0; j < n_bas_fcts; j++)
1190 {
1191 phi_b = PHI(bas_fcts, j, bm[i]); /* !!! phi[j](bm[i]); */
1192 for (k=0; k<DIM_OF_WORLD; k++)
1193 vd[k] += el_vec_d[j][k]*phi_b;
1194 }
1195 vm[i] = NORM_DOW(vd);
1196 }
1197
1198 if (el_parametric)
1199 {
1200 el_parametric->coord_to_world(el_info, NULL, N_EDGES_2D,
1201 (const REAL_B *) bm, xm);
1202 }
1203 else
1204 {
1205 for (j = 0; j < DIM_OF_WORLD; j++)
1206 {
1207 xm[0][j] = 0.5*(x[1][j] + x[2][j]);
1208 xm[1][j] = 0.5*(x[0][j] + x[2][j]);
1209 xm[2][j] = 0.5*(x[0][j] + x[1][j]);
1210 }
1211 }
1212
1213 bnew[0] = b[0]; bnew[1] = bm[2]; bnew[2] = bm[1];
1214 xnew[0] = x[0]; xnew[1] = xm[2]; xnew[2] = xm[1];
1215 vnew[0] = v[0]; vnew[1] = vm[2]; vnew[2] = vm[1];
1216 graph_value_d_recursive(refine-1, bnew, vnew, xnew);
1217
1218 bnew[0] = bm[2]; bnew[1] = b[1]; bnew[2] = bm[0];
1219 xnew[0] = xm[2]; xnew[1] = x[1]; xnew[2] = xm[0];
1220 vnew[0] = vm[2]; vnew[1] = v[1]; vnew[2] = vm[0];
1221 graph_value_d_recursive(refine-1, bnew, vnew, xnew);
1222
1223 bnew[0] = b[2]; bnew[1] = bm[1]; bnew[2] = bm[0];
1224 xnew[0] = x[2]; xnew[1] = xm[1]; xnew[2] = xm[0];
1225 vnew[0] = v[2]; vnew[1] = vm[1]; vnew[2] = vm[0];
1226 graph_value_d_recursive(refine-1, bnew, vnew, xnew);
1227
1228 bnew[0] = bm[0]; bnew[1] = bm[1]; bnew[2] = bm[2];
1229 xnew[0] = xm[0]; xnew[1] = xm[1]; xnew[2] = xm[2];
1230 vnew[0] = vm[0]; vnew[1] = vm[1]; vnew[2] = vm[2];
1231 graph_value_d_recursive(refine-1, bnew, vnew, xnew);
1232
1233 return;
1234 }
1235
1236 /* draw triangle */
1237
1238 glBegin(GL_TRIANGLES);
1239 for (i = 0; i < N_VERTICES_2D; i++) {
1240 glColor3fv(val_color(v[i]));
1241 glVertex2dv(x[i]);
1242 }
1243
1244 glEnd();
1245 }
1246
1247
graph_value_d_fct(const EL_INFO * elinfo,void * data)1248 static void graph_value_d_fct(const EL_INFO *elinfo, void *data)
1249 {
1250 FUNCNAME("graph_value_d_fct");
1251 static REAL_B bnew[3]={INIT_BARY_2D(1.0,0.0,0.0),
1252 INIT_BARY_2D(0.0,1.0,0.0),
1253 INIT_BARY_2D(0.0,0.0,1.0)};
1254 static const REAL *b[3]={bnew[0],bnew[1],bnew[2]};
1255 static const REAL *x[3];
1256 static REAL_D coord[N_VERTICES_2D];
1257 static REAL v[3], phi_b;
1258 static REAL_D vd;
1259 int i, j, k;
1260 PARAMETRIC *parametric = elinfo->mesh->parametric;
1261
1262 if (bas_fcts->get_real_d_vec)
1263 el_vec_d = fill_el_real_d_vec(NULL, elinfo->el, drdv)->vec;
1264 else
1265 ERROR("no bas_fcts->get_real_d_vec()\n");
1266
1267 el_info = elinfo;
1268 if (parametric)
1269 {
1270 el_parametric = parametric;
1271 parametric->init_element(elinfo, parametric);
1272 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D, vertex_lambda_2d,
1273 coord);
1274 for (i = 0; i < 3; i++)
1275 x[i] = coord[i];
1276 }
1277 else
1278 {
1279 el_parametric = NULL;
1280 for (i = 0; i < 3; i++)
1281 x[i] = elinfo->coord[i];
1282 }
1283
1284 for (i = 0; i < 3; i++)
1285 {
1286 SET_DOW(0.0, vd);
1287 for (j = 0; j < n_bas_fcts; j++)
1288 {
1289 phi_b = PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
1290 for (k = 0; k < DIM_OF_WORLD; k++)
1291 vd[k] += el_vec_d[j][k]*phi_b;
1292 }
1293 v[i] = NORM_DOW(vd);
1294 }
1295
1296 graph_value_d_recursive(nrefine, b, v, x);
1297 }
1298
1299
graph_drv_d_2d(GRAPH_WINDOW win,const DOF_REAL_D_VEC * v,REAL min,REAL max,int refine)1300 static void graph_drv_d_2d(GRAPH_WINDOW win, const DOF_REAL_D_VEC *v,
1301 REAL min, REAL max, int refine)
1302 {
1303 FUNCNAME("graph_drv_d_2d");
1304 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
1305
1306 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
1307 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
1308 return;
1309 }
1310
1311 drdv = v;
1312 bas_fcts = v->fe_space->bas_fcts;
1313 n_bas_fcts = bas_fcts->n_bas_fcts;
1314 phi = bas_fcts->phi;
1315
1316 if (refine >= 0)
1317 nrefine = refine;
1318 else
1319 nrefine = MAX(0, bas_fcts->degree-1);
1320
1321 if (min < max) {
1322 val_min = min;
1323 val_max = max;
1324 }
1325 else {
1326 val_max = -(val_min = 1.0E20);
1327 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL,
1328 val_d_minmax_fct, NULL);
1329 val_max = MAX(val_max, val_min + 1.0E-10);
1330 }
1331 val_fac = 1.0 / (val_max - val_min);
1332
1333
1334 OGL_set_std_window(ogl_win);
1335
1336 mesh_traverse(v->fe_space->admin->mesh, -1, CALL_LEAF_EL|FILL_COORDS,
1337 graph_value_d_fct, NULL);
1338
1339 OGL_FLUSH(ogl_win);
1340 }
1341
1342
1343 /****************************************************************************/
1344 /* graph_fvalue: general graph_value without DOF_REAL_VEC */
1345 /****************************************************************************/
1346 static REAL (*fvalue_fct)(const EL_INFO *, const REAL *lambda) = NULL;
1347
gval_minmax_fct(const EL_INFO * elinfo,void * data)1348 static void gval_minmax_fct(const EL_INFO *elinfo, void *data)
1349 {
1350 FUNCNAME("gval_minmax_fct");
1351 #define NTEST 4
1352 static REAL_B b[NTEST] = {INIT_BARY_2D(1.0, 0.0, 0.0),
1353 INIT_BARY_2D(0.0, 1.0, 0.0),
1354 INIT_BARY_2D(0.0, 0.0, 1.0),
1355 INIT_BARY_2D(0.333,0.333,0.334)};
1356 int i, ntest;
1357 REAL v;
1358
1359 TEST_EXIT(fvalue_fct, "no fvalue_fct\n");
1360
1361 if (nrefine>0) ntest = NTEST;
1362 else ntest = 3;
1363
1364 for (i = 0; i < ntest; i++)
1365 {
1366 v = fvalue_fct(elinfo, b[i]);
1367
1368 val_min = MIN(val_min, v);
1369 val_max = MAX(val_max, v);
1370 }
1371 #undef NTEST
1372 }
1373
1374 /****************************************************************************/
1375
graph_fvalue_recursive(int refine,const REAL * b[3],REAL v[3],const REAL * x[3])1376 static void graph_fvalue_recursive(int refine, const REAL *b[3],
1377 REAL v[3], const REAL *x[3])
1378 {
1379 int i,j;
1380
1381 if (refine > 0) /* refine and call recursively */
1382 {
1383 const REAL *bnew[3], *xnew[3];
1384 REAL_B bm[3];
1385 REAL vnew[3], vm[3], xm[3][DIM_OF_WORLD];
1386
1387 for (j = 0; j < 3; j++)
1388 {
1389 bm[0][j] = 0.5*(b[1][j] + b[2][j]);
1390 bm[1][j] = 0.5*(b[0][j] + b[2][j]);
1391 bm[2][j] = 0.5*(b[0][j] + b[1][j]);
1392 }
1393
1394 for (i = 0; i < 3; i++)
1395 {
1396 vm[i] = fvalue_fct(el_info, bm[i]);
1397 }
1398
1399 if (el_parametric)
1400 {
1401 el_parametric->coord_to_world(el_info, NULL, N_EDGES_2D,
1402 (const REAL_B *) bm, xm);
1403 }
1404 else
1405 {
1406 for (j = 0; j < DIM_OF_WORLD; j++)
1407 {
1408 xm[0][j] = 0.5*(x[1][j] + x[2][j]);
1409 xm[1][j] = 0.5*(x[0][j] + x[2][j]);
1410 xm[2][j] = 0.5*(x[0][j] + x[1][j]);
1411 }
1412 }
1413
1414 bnew[0] = b[0]; bnew[1] = bm[2]; bnew[2] = bm[1];
1415 xnew[0] = x[0]; xnew[1] = xm[2]; xnew[2] = xm[1];
1416 vnew[0] = v[0]; vnew[1] = vm[2]; vnew[2] = vm[1];
1417 graph_fvalue_recursive(refine-1, bnew, vnew, xnew);
1418
1419 bnew[0] = bm[2]; bnew[1] = b[2]; bnew[2] = bm[1];
1420 xnew[0] = xm[2]; xnew[1] = x[2]; xnew[2] = xm[1];
1421 vnew[0] = vm[2]; vnew[1] = v[2]; vnew[2] = vm[1];
1422 graph_fvalue_recursive(refine-1, bnew, vnew, xnew);
1423
1424 bnew[0] = b[2]; bnew[1] = bm[2]; bnew[2] = bm[0];
1425 xnew[0] = x[2]; xnew[1] = xm[2]; xnew[2] = xm[0];
1426 vnew[0] = v[2]; vnew[1] = vm[2]; vnew[2] = vm[0];
1427 graph_fvalue_recursive(refine-1, bnew, vnew, xnew);
1428
1429 bnew[0] = bm[2]; bnew[1] = b[1]; bnew[2] = bm[0];
1430 xnew[0] = xm[2]; xnew[1] = x[1]; xnew[2] = xm[0];
1431 vnew[0] = vm[2]; vnew[1] = v[1]; vnew[2] = vm[0];
1432 graph_fvalue_recursive(refine-1, bnew, vnew, xnew);
1433
1434 return;
1435 }
1436
1437 /* draw triangle */
1438
1439 glBegin(GL_TRIANGLES);
1440 for (i = 0; i < N_VERTICES_2D; i++) {
1441 glColor3fv(val_color(v[i]));
1442 glVertex2dv(x[i]);
1443 }
1444
1445 glEnd();
1446 }
1447
1448
graph_fvalue_fct(const EL_INFO * elinfo,void * data)1449 static void graph_fvalue_fct(const EL_INFO *elinfo, void *data)
1450 {
1451 FUNCNAME("graph_fvalue_fct");
1452 static REAL_B bnew[3]={INIT_BARY_2D(1.0,0.0,0.0),
1453 INIT_BARY_2D(0.0,1.0,0.0),
1454 INIT_BARY_2D(0.0,0.0,1.0)};
1455 static const REAL *b[3]={bnew[0],bnew[1],bnew[2]};
1456 static const REAL *x[3];
1457 static REAL_D coord[N_VERTICES_2D];
1458 static REAL v[3];
1459 int i;
1460 PARAMETRIC *parametric = elinfo->mesh->parametric;
1461
1462 TEST_EXIT(fvalue_fct, "no fvalue_fct\n");
1463
1464 el_info = elinfo;
1465 if (parametric)
1466 {
1467 el_parametric = parametric;
1468 parametric->init_element(elinfo, parametric);
1469 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D, vertex_lambda_2d,
1470 coord);
1471 for (i = 0; i < 3; i++)
1472 x[i] = coord[i];
1473 }
1474 else
1475 {
1476 el_parametric = NULL;
1477 for (i = 0; i < 3; i++)
1478 x[i] = elinfo->coord[i];
1479 }
1480
1481 for (i = 0; i < 3; i++)
1482 {
1483 v[i] = fvalue_fct(elinfo, b[i]);
1484 }
1485
1486 graph_fvalue_recursive(nrefine, b, v, x);
1487 }
1488
graph_fvalues_2d(GRAPH_WINDOW win,MESH * mesh,REAL (* fct)(const EL_INFO *,const REAL * lambda),FLAGS fill_flag,REAL min,REAL max,int refine)1489 void graph_fvalues_2d(GRAPH_WINDOW win, MESH *mesh,
1490 REAL(*fct)(const EL_INFO *, const REAL *lambda),
1491 FLAGS fill_flag, REAL min, REAL max, int refine)
1492 {
1493 FUNCNAME("graph_fvalues_2d");
1494 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
1495
1496 TEST_EXIT(mesh, "no mesh\n");
1497 TEST_EXIT(fvalue_fct = fct, "no fct\n");
1498
1499 if((DIM_OF_WORLD != 2) || (mesh->dim != 2)) {
1500 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
1501 return;
1502 }
1503
1504 nrefine = MAX(refine, 0);
1505
1506 if (min < max) {
1507 val_min = min;
1508 val_max = max;
1509 }
1510 else {
1511 val_max = -(val_min = 1.0E20);
1512 mesh_traverse(mesh, -1, CALL_LEAF_EL | fill_flag, gval_minmax_fct, NULL);
1513 val_max = MAX(val_max, val_min + 1.0E-10);
1514 }
1515 val_fac = 1.0 / (val_max - val_min);
1516
1517 OGL_set_std_window(ogl_win);
1518
1519 mesh_traverse(mesh, -1, CALL_LEAF_EL|FILL_COORDS, graph_fvalue_fct, NULL);
1520
1521 OGL_FLUSH(ogl_win);
1522 }
1523
1524 /****************************************************************************/
1525 /****************************************************************************/
1526 /* MULTIGRID LEVEL DISPLAYS */
1527 /****************************************************************************/
1528 /****************************************************************************/
1529
graph_mesh_mg_2d(GRAPH_WINDOW win,MESH * mesh,const GRAPH_RGBCOLOR c,FLAGS flags,int mg_level)1530 void graph_mesh_mg_2d(GRAPH_WINDOW win, MESH *mesh, const GRAPH_RGBCOLOR c,
1531 FLAGS flags, int mg_level)
1532 {
1533 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
1534
1535 if (!mesh) return;
1536
1537 if((DIM_OF_WORLD != 2) || (mesh->dim != 2)) {
1538 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
1539 return;
1540 }
1541
1542 OGL_set_std_window(ogl_win);
1543 if (c) linecolor_2d = c;
1544 else linecolor_2d = rgb_black;
1545 graph_mesh_flags_2d = flags;
1546 mesh_traverse(mesh, mg_level, CALL_MG_LEVEL | FILL_COORDS | FILL_BOUND,
1547 OGLgraph_mesh_fct_2d, NULL);
1548
1549 OGL_FLUSH(ogl_win);
1550
1551 return;
1552 }
1553
1554 /****************************************************************************/
1555 static const int *mg_sort_dof_invers = NULL;
1556 static int mg_el_vec_size = 0;
1557 static REAL *mg_el_vec = NULL;
1558 static const REAL *drv_vec = NULL;
1559 static const DOF_ADMIN *mg_admin = NULL;
1560
1561
val_minmax_mg_fct(const EL_INFO * elinfo,void * data)1562 static void val_minmax_mg_fct(const EL_INFO *elinfo, void *data)
1563 {
1564 FUNCNAME("val_minmax_mg_fct");
1565 #define NTEST 4
1566 static REAL_B b[NTEST] = {INIT_BARY_2D(1.0, 0.0, 0.0),
1567 INIT_BARY_2D(0.0, 1.0, 0.0),
1568 INIT_BARY_2D(0.0, 0.0, 1.0),
1569 INIT_BARY_2D(0.333,0.333,0.334)};
1570 int i, j, ntest;
1571 REAL v;
1572 const DOF *dof_indices;
1573
1574 if (bas_fcts->get_dof_indices)
1575 dof_indices = GET_DOF_INDICES(bas_fcts, elinfo->el, mg_admin, NULL)->vec;
1576 else {
1577 ERROR("no bas_fcts->get_dof_indices()\n");
1578 return;
1579 }
1580
1581 if (mg_sort_dof_invers) {
1582 for (i = 0; i < n_bas_fcts; i++) {
1583 mg_el_vec[i] = drv_vec[mg_sort_dof_invers[dof_indices[i]]];
1584 }
1585 }
1586 else /* no mg_sort_dof_invers */
1587 {
1588 for (i = 0; i < n_bas_fcts; i++) {
1589 mg_el_vec[i] = drv_vec[dof_indices[i]];
1590 }
1591 }
1592
1593 if (nrefine>0) ntest = NTEST;
1594 else ntest = 3;
1595
1596 for (i = 0; i < ntest; i++)
1597 {
1598 v = 0.0;
1599 for (j=0; j<n_bas_fcts; j++)
1600 v += mg_el_vec[j] * PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
1601
1602 val_min = MIN(val_min, v);
1603 val_max = MAX(val_max, v);
1604 }
1605 #undef NTEST
1606 }
1607
graph_value_mg_recursive(int refine,const REAL * b[3],REAL v[3],const REAL * x[3])1608 static void graph_value_mg_recursive(int refine, const REAL *b[3],
1609 REAL v[3], const REAL *x[3])
1610 {
1611 int i, j;
1612
1613 if (refine > 0) /* refine and call recursively */
1614 {
1615 const REAL *bnew[3], *xnew[3];
1616 REAL_B bm[3];
1617 REAL vnew[3], vm[3], xm[3][DIM_OF_WORLD];
1618
1619 for (j = 0; j < 3; j++)
1620 {
1621 bm[0][j] = 0.5*(b[1][j] + b[2][j]);
1622 bm[1][j] = 0.5*(b[0][j] + b[2][j]);
1623 bm[2][j] = 0.5*(b[0][j] + b[1][j]);
1624 }
1625
1626 for (i = 0; i < 3; i++)
1627 {
1628 vm[i] = 0.0;
1629 for (j=0; j<n_bas_fcts; j++)
1630 vm[i] += mg_el_vec[j] * PHI(bas_fcts, j, bm[i]); /* !!! phi[j](bm[i]); */
1631 }
1632
1633 if (el_parametric)
1634 {
1635 el_parametric->coord_to_world(el_info, NULL, N_EDGES_2D,
1636 (const REAL_B *) bm, xm);
1637 }
1638 else
1639 {
1640 for (j = 0; j < DIM_OF_WORLD; j++)
1641 {
1642 xm[0][j] = 0.5*(x[1][j] + x[2][j]);
1643 xm[1][j] = 0.5*(x[0][j] + x[2][j]);
1644 xm[2][j] = 0.5*(x[0][j] + x[1][j]);
1645 }
1646 }
1647
1648 bnew[0] = b[0]; bnew[1] = bm[2]; bnew[2] = bm[1];
1649 xnew[0] = x[0]; xnew[1] = xm[2]; xnew[2] = xm[1];
1650 vnew[0] = v[0]; vnew[1] = vm[2]; vnew[2] = vm[1];
1651 graph_value_mg_recursive(refine-1, bnew, vnew, xnew);
1652
1653 bnew[0] = bm[2]; bnew[1] = b[2]; bnew[2] = bm[1];
1654 xnew[0] = xm[2]; xnew[1] = x[2]; xnew[2] = xm[1];
1655 vnew[0] = vm[2]; vnew[1] = v[2]; vnew[2] = vm[1];
1656 graph_value_mg_recursive(refine-1, bnew, vnew, xnew);
1657
1658 bnew[0] = b[2]; bnew[1] = bm[2]; bnew[2] = bm[0];
1659 xnew[0] = x[2]; xnew[1] = xm[2]; xnew[2] = xm[0];
1660 vnew[0] = v[2]; vnew[1] = vm[2]; vnew[2] = vm[0];
1661 graph_value_mg_recursive(refine-1, bnew, vnew, xnew);
1662
1663 bnew[0] = bm[2]; bnew[1] = b[1]; bnew[2] = bm[0];
1664 xnew[0] = xm[2]; xnew[1] = x[1]; xnew[2] = xm[0];
1665 vnew[0] = vm[2]; vnew[1] = v[1]; vnew[2] = vm[0];
1666 graph_value_mg_recursive(refine-1, bnew, vnew, xnew);
1667
1668 return;
1669 }
1670
1671 /* draw triangle */
1672
1673 glBegin(GL_TRIANGLES);
1674 for (i = 0; i < N_VERTICES_2D; i++) {
1675 glColor3fv(val_color(v[i]));
1676 glVertex2dv(x[i]);
1677 }
1678
1679 glEnd();
1680 }
1681
1682
graph_value_mg_fct(const EL_INFO * elinfo,void * data)1683 static void graph_value_mg_fct(const EL_INFO *elinfo, void *data)
1684 {
1685 FUNCNAME("graph_value_mg_fct");
1686 static REAL_B bnew[3]={INIT_BARY_2D(1.0,0.0,0.0),
1687 INIT_BARY_2D(0.0,1.0,0.0),
1688 INIT_BARY_2D(0.0,0.0,1.0)};
1689 static const REAL *b[3]={bnew[0],bnew[1],bnew[2]};
1690 const REAL *x[3];
1691 REAL_D coord[N_VERTICES_2D];
1692 REAL v[3];
1693 int i, j;
1694 const DOF *dof_indices;
1695 PARAMETRIC *parametric = elinfo->mesh->parametric;
1696
1697 if (bas_fcts->get_dof_indices)
1698 dof_indices = GET_DOF_INDICES(bas_fcts, elinfo->el, mg_admin, NULL)->vec;
1699 else {
1700 ERROR("no bas_fcts->get_dof_indices()\n");
1701 return;
1702 }
1703
1704 if (mg_sort_dof_invers) {
1705 for (i = 0; i < n_bas_fcts; i++) {
1706 mg_el_vec[i] = drv_vec[mg_sort_dof_invers[dof_indices[i]]];
1707 }
1708 }
1709 else /* no mg_sort_dof_invers */
1710 {
1711 for (i = 0; i < n_bas_fcts; i++) {
1712 mg_el_vec[i] = drv_vec[dof_indices[i]];
1713 }
1714 }
1715
1716 el_info = elinfo;
1717 if (parametric)
1718 {
1719 el_parametric = parametric;
1720 parametric->init_element(elinfo, parametric);
1721 parametric->coord_to_world(elinfo, NULL, N_VERTICES_2D, vertex_lambda_2d,
1722 coord);
1723 for (i = 0; i < 3; i++)
1724 x[i] = coord[i];
1725 }
1726 else
1727 {
1728 el_parametric = NULL;
1729 for (i = 0; i < 3; i++)
1730 x[i] = elinfo->coord[i];
1731 }
1732
1733 for (i = 0; i < 3; i++)
1734 {
1735 v[i] = 0.0;
1736 for (j=0; j<n_bas_fcts; j++)
1737 v[i] += mg_el_vec[j] * PHI(bas_fcts, j, b[i]); /* !!! phi[j](b[i]); */
1738 }
1739
1740 graph_value_mg_recursive(nrefine, b, v, x);
1741 }
1742
1743
graph_values_mg_2d(GRAPH_WINDOW win,const DOF_REAL_VEC * v,REAL min,REAL max,int refine,int mg_level,const FE_SPACE * fe_space,const int * sort_dof_invers)1744 void graph_values_mg_2d(GRAPH_WINDOW win, const DOF_REAL_VEC *v,
1745 REAL min, REAL max, int refine,
1746 int mg_level, const FE_SPACE *fe_space,
1747 const int *sort_dof_invers)
1748 {
1749 FUNCNAME("graph_values_mg_2d");
1750 OGL_WINDOW *ogl_win = (OGL_WINDOW *)win;
1751
1752 TEST_EXIT(v && fe_space && fe_space->admin,
1753 "no vec or fe_space or admin\n");
1754
1755 if((DIM_OF_WORLD != 2) || (v->fe_space->admin->mesh->dim != 2)) {
1756 ERROR("Only implemented for DIM_OF_WORLD==2 and dim==2!\n");
1757 return;
1758 }
1759
1760 drv = v;
1761 drv_vec = drv->vec;
1762 bas_fcts = fe_space->bas_fcts;
1763 mg_admin = fe_space->admin;
1764 n_bas_fcts = bas_fcts->n_bas_fcts;
1765 phi = bas_fcts->phi;
1766
1767 mg_sort_dof_invers = sort_dof_invers;
1768 if (n_bas_fcts > mg_el_vec_size) {
1769 mg_el_vec = MEM_REALLOC(mg_el_vec, mg_el_vec_size, n_bas_fcts, REAL);
1770 mg_el_vec_size = n_bas_fcts;
1771 }
1772
1773 if (refine >= 0)
1774 nrefine = refine;
1775 else
1776 nrefine = MAX(0, bas_fcts->degree-1);
1777
1778 if (min < max) {
1779 val_min = min;
1780 val_max = max;
1781 }
1782 else {
1783 val_max = -(val_min = 1.0E20);
1784 mesh_traverse(fe_space->admin->mesh, mg_level,
1785 CALL_MG_LEVEL, val_minmax_mg_fct, NULL);
1786 MSG("<%s> value range in [%.3le , %.3le]\n", NAME(drv), val_min, val_max);
1787
1788 val_max = MAX(val_max, val_min + 1.0E-5);
1789 }
1790 val_fac = 1.0 / (val_max - val_min);
1791
1792
1793 OGL_set_std_window(ogl_win);
1794
1795 mesh_traverse(fe_space->admin->mesh, mg_level, CALL_MG_LEVEL | FILL_COORDS,
1796 graph_value_mg_fct, NULL);
1797
1798 graph_mesh_mg_2d(win, fe_space->admin->mesh, rgb_black, 0, mg_level);
1799
1800 OGL_FLUSH(ogl_win);
1801 }
1802