Codebase list gerris / 8731f8d
MetricCubed takes an optional "level" parameter darcs-hash:20091028001441-d4795-9ece460adb33c506f2c5d3cec027b0ab2f1c1a35.gz Stephane Popinet 14 years ago
4 changed file(s) with 144 addition(s) and 42 deletion(s). Raw diff Collapse all Expand all
+3
-0
src/domain.c less more
26612661 g_return_if_fail (child.c[n]->data == NULL);
26622662 child.c[n]->data = g_malloc0 (gfs_domain_variables_size (domain));
26632663 }
2664 if (GFS_CELL_IS_BOUNDARY (cell))
2665 for (n = 0; n < FTT_CELLS; n++)
2666 child.c[n]->flags |= GFS_FLAG_BOUNDARY;
26642667 }
26652668 }
26662669
+137
-42
src/metric.c less more
1616 * 02111-1307, USA.
1717 */
1818
19 #include <stdlib.h>
1920 #include "metric.h"
2021 #include <complex.h>
2122 #include "map.h"
368369 {
369370 }
370371
372 typedef struct {
373 double x, y, x1, y1, z1;
374 double a;
375 } Point;
376
377 typedef struct {
378 Point * v[4];
379 double a, h[4];
380 } Square;
381
382 static void point_new (Point * p, double x, double y)
383 {
384 p->x = x; p->y = y;
385 cmap_xy2XYZ (x, y, &p->x1, &p->y1, &p->z1);
386 }
387
388 static Point ** matrix_refine (Point ** m, int n)
389 {
390 int n1 = 2*n - 1, i, j;
391 Point ** r = gfs_matrix_new (n1, n1, sizeof (Point));
392 for (i = 0; i < n; i++)
393 for (j = 0; j < n; j++)
394 r[2*i][2*j] = m[i][j];
395 for (i = 0; i < n - 1; i++)
396 for (j = 0; j < n - 1; j++) {
397 point_new (&r[2*i+1][2*j], (m[i][j].x + m[i+1][j].x)/2., m[i][j].y);
398 point_new (&r[2*i][2*j+1], m[i][j].x, (m[i][j].y + m[i][j+1].y)/2.);
399 point_new (&r[2*i+1][2*j+1], (m[i][j].x + m[i+1][j].x)/2., (m[i][j].y + m[i][j+1].y)/2.);
400 }
401 i = n - 1;
402 for (j = 0; j < n - 1; j++)
403 point_new (&r[2*i][2*j+1], m[i][j].x, (m[i][j].y + m[i][j+1].y)/2.);
404 j = n - 1;
405 for (i = 0; i < n - 1; i++)
406 point_new (&r[2*i+1][2*j], (m[i][j].x + m[i+1][j].x)/2., m[i][j].y);
407 gfs_matrix_free (m);
408 return r;
409 }
410
411 static Point ** matrix_from_cell (FttCell * cell)
412 {
413 FttVector p;
414 ftt_cell_pos (cell, &p);
415 double h = ftt_cell_size (cell)/2.;
416 Point ** r = gfs_matrix_new (2, 2, sizeof (Point));
417 point_new (&r[0][0], p.x - h, p.y - h);
418 point_new (&r[1][0], p.x + h, p.y - h);
419 point_new (&r[1][1], p.x + h, p.y + h);
420 point_new (&r[0][1], p.x - h, p.y + h);
421 return r;
422 }
423
424 static double matrix_a (Point ** r, int m, int i0, int j0)
425 {
426 int i, j;
427 double a = 0.;
428 double h = r[m][0].x - r[0][0].x;
429 for (i = 0; i < m; i++)
430 for (j = 0; j < m; j++)
431 a += excess_of_quad (&r[i0+i][j0+j].x1, &r[i0+i+1][j0+j].x1,
432 &r[i0+i+1][j0+j+1].x1, &r[i0+i][j0+j+1].x1);
433 return 4.*a/(M_PI*M_PI*h*h);
434 }
435
436 static double matrix_hx (Point ** r, int m, int i0, int j0)
437 {
438 int i;
439 double hx = 0.;
440 double h = r[m][0].x - r[0][0].x;
441 for (i = 0; i < m; i++)
442 hx += angle_between_vectors (&r[i0+i][j0].x1, &r[i0+i+1][j0].x1);
443 return 2.*hx/(M_PI*h);
444 }
445
446 static double matrix_hy (Point ** r, int m, int i0, int j0)
447 {
448 int j;
449 double hy = 0.;
450 double h = r[m][0].x - r[0][0].x;
451 for (j = 0; j < m; j++)
452 hy += angle_between_vectors (&r[i0][j0+j].x1, &r[i0][j0+j+1].x1);
453 return 2.*hy/(M_PI*h);
454 }
455
371456 static void cubed_coarse_fine (FttCell * parent, GfsVariable * a)
372457 {
458 if (GFS_CELL_IS_BOUNDARY (parent))
459 return;
460
373461 GfsMetricCubed * cubed = GTS_OBJECT (a)->reserved;
462 g_assert (GFS_IS_METRIC_CUBED (cubed));
463 Point ** r = matrix_from_cell (parent);
464 r = matrix_refine (r, 2);
465 int n = 3, level = cubed->level - (ftt_cell_level (parent) + 1);
466 while (level-- > 0) {
467 r = matrix_refine (r, n);
468 n = 2*n - 1;
469 }
470
374471 FttCellChildren child;
375472 ftt_cell_children (parent, &child);
376 FttVector p;
377 ftt_cell_pos (parent, &p);
378 double h = ftt_cell_size (parent)/2.;
379 double v0[3];
380 v0[0] = p.x; v0[1] = p.y;
381 cmap_xy2XYZ (v0[0], v0[1], &v0[0], &v0[1], &v0[2]);
382 double v1[3], v2[3], v3[3], v4[3];
383 v1[0] = p.x - h; v1[1] = p.y - h;
384 cmap_xy2XYZ (v1[0], v1[1], &v1[0], &v1[1], &v1[2]);
385 v2[0] = p.x - h; v2[1] = p.y + h;
386 cmap_xy2XYZ (v2[0], v2[1], &v2[0], &v2[1], &v2[2]);
387 v3[0] = p.x + h; v3[1] = p.y + h;
388 cmap_xy2XYZ (v3[0], v3[1], &v3[0], &v3[1], &v3[2]);
389 v4[0] = p.x + h; v4[1] = p.y - h;
390 cmap_xy2XYZ (v4[0], v4[1], &v4[0], &v4[1], &v4[2]);
391 double v5[3], v6[3], v7[3], v8[3];
392 v5[0] = p.x; v5[1] = p.y - h;
393 cmap_xy2XYZ (v5[0], v5[1], &v5[0], &v5[1], &v5[2]);
394 v6[0] = p.x - h; v6[1] = p.y;
395 cmap_xy2XYZ (v6[0], v6[1], &v6[0], &v6[1], &v6[2]);
396 v7[0] = p.x; v7[1] = p.y + h;
397 cmap_xy2XYZ (v7[0], v7[1], &v7[0], &v7[1], &v7[2]);
398 v8[0] = p.x + h; v8[1] = p.y;
399 cmap_xy2XYZ (v8[0], v8[1], &v8[0], &v8[1], &v8[2]);
400
401 GFS_VALUE (child.c[0], a) = 4.*excess_of_quad (v6, v0, v7, v2)/(M_PI*M_PI*h*h);
402 GFS_VALUE (child.c[1], a) = 4.*excess_of_quad (v0, v8, v3, v7)/(M_PI*M_PI*h*h);
403 GFS_VALUE (child.c[2], a) = 4.*excess_of_quad (v1, v5, v0, v6)/(M_PI*M_PI*h*h);
404 GFS_VALUE (child.c[3], a) = 4.*excess_of_quad (v5, v4, v8, v0)/(M_PI*M_PI*h*h);
473 int m = n/2;
474
475 GFS_VALUE (child.c[0], a) = matrix_a (r, m, 0, m);
476 GFS_VALUE (child.c[1], a) = matrix_a (r, m, m, m);
477 GFS_VALUE (child.c[2], a) = matrix_a (r, m, 0, 0);
478 GFS_VALUE (child.c[3], a) = matrix_a (r, m, m, 0);
405479
406480 GFS_VALUE (child.c[0], cubed->h[0]) = GFS_VALUE (child.c[1], cubed->h[1]) =
407 2.*angle_between_vectors (v0, v7)/(M_PI*h);
481 matrix_hy (r, m, m, m);
408482 GFS_VALUE (child.c[0], cubed->h[3]) = GFS_VALUE (child.c[2], cubed->h[2]) =
409 2.*angle_between_vectors (v0, v6)/(M_PI*h);
483 matrix_hx (r, m, 0, m);
410484 GFS_VALUE (child.c[2], cubed->h[0]) = GFS_VALUE (child.c[3], cubed->h[1]) =
411 2.*angle_between_vectors (v0, v5)/(M_PI*h);
485 matrix_hy (r, m, m, 0);
412486 GFS_VALUE (child.c[1], cubed->h[3]) = GFS_VALUE (child.c[3], cubed->h[2]) =
413 2.*angle_between_vectors (v0, v8)/(M_PI*h);
414
415 GFS_VALUE (child.c[0], cubed->h[2]) = 2.*angle_between_vectors (v2, v7)/(M_PI*h);
416 GFS_VALUE (child.c[0], cubed->h[1]) = 2.*angle_between_vectors (v2, v6)/(M_PI*h);
417 GFS_VALUE (child.c[1], cubed->h[2]) = 2.*angle_between_vectors (v3, v7)/(M_PI*h);
418 GFS_VALUE (child.c[1], cubed->h[0]) = 2.*angle_between_vectors (v3, v8)/(M_PI*h);
419 GFS_VALUE (child.c[2], cubed->h[3]) = 2.*angle_between_vectors (v1, v5)/(M_PI*h);
420 GFS_VALUE (child.c[2], cubed->h[1]) = 2.*angle_between_vectors (v1, v6)/(M_PI*h);
421 GFS_VALUE (child.c[3], cubed->h[0]) = 2.*angle_between_vectors (v4, v8)/(M_PI*h);
422 GFS_VALUE (child.c[3], cubed->h[3]) = 2.*angle_between_vectors (v4, v5)/(M_PI*h);
487 matrix_hx (r, m, m, m);
488
489 GFS_VALUE (child.c[0], cubed->h[2]) = matrix_hx (r, m, 0, n - 1);
490 GFS_VALUE (child.c[0], cubed->h[1]) = matrix_hy (r, m, 0, m);
491 GFS_VALUE (child.c[1], cubed->h[2]) = matrix_hx (r, m, m, n - 1);
492 GFS_VALUE (child.c[1], cubed->h[0]) = matrix_hy (r, m, n - 1, m);
493 GFS_VALUE (child.c[2], cubed->h[3]) = matrix_hx (r, m, 0, 0);
494 GFS_VALUE (child.c[2], cubed->h[1]) = matrix_hy (r, m, 0, 0);
495 GFS_VALUE (child.c[3], cubed->h[0]) = matrix_hy (r, m, n - 1, 0);
496 GFS_VALUE (child.c[3], cubed->h[3]) = matrix_hx (r, m, m, 0);
497
498 gfs_matrix_free (r);
423499 }
424500
425501 static void cubed_fine_coarse (FttCell * parent, GfsVariable * a)
426502 {
427503 GfsMetricCubed * cubed = GTS_OBJECT (a)->reserved;
504 g_assert (GFS_IS_METRIC_CUBED (cubed));
428505 FttCellChildren child;
429506 guint n;
430507
444521 GFS_VALUE (child.c[3], cubed->h[3]))/2.;
445522 }
446523
524 static void metric_cubed_write (GtsObject * o, FILE * fp)
525 {
526 (* GTS_OBJECT_CLASS (gfs_metric_cubed_class ())->parent_class->write) (o, fp);
527 if (GFS_METRIC_CUBED (o)->level != 0)
528 fprintf (fp, " %d", GFS_METRIC_CUBED (o)->level);
529 }
530
447531 static void metric_cubed_read (GtsObject ** o, GtsFile * fp)
448532 {
449533 (* GTS_OBJECT_CLASS (gfs_metric_cubed_class ())->parent_class->read) (o, fp);
457541 }
458542
459543 GfsMetricCubed * cubed = GFS_METRIC_CUBED (*o);
544 if (fp->type == GTS_INT) {
545 cubed->level = atoi (fp->token->str);
546 gts_file_next_token (fp);
547 }
548
460549 FttDirection d;
461550 for (d = 0; d < FTT_NEIGHBORS; d++) {
462551 gchar * name = g_strdup_printf ("Ch%d", d);
482571 static void metric_cubed_class_init (GtsObjectClass * klass)
483572 {
484573 klass->read = metric_cubed_read;
574 klass->write = metric_cubed_write;
485575 }
486576
487577 static void metric_cubed_init (GfsEvent * m)
623713
624714 static void lonlat_coarse_fine (FttCell * parent, GfsVariable * a)
625715 {
716 if (GFS_CELL_IS_BOUNDARY (parent))
717 return;
718
626719 GfsMetricLonLat * lonlat = GTS_OBJECT (a)->reserved;
720 g_assert (GFS_IS_METRIC_LON_LAT (lonlat));
627721 FttCellChildren child;
628722 ftt_cell_children (parent, &child);
629723 FttVector p;
649743 static void lonlat_fine_coarse (FttCell * parent, GfsVariable * a)
650744 {
651745 GfsMetricLonLat * lonlat = GTS_OBJECT (a)->reserved;
746 g_assert (GFS_IS_METRIC_LON_LAT (lonlat));
652747 FttCellChildren child;
653748 guint n;
654749
+1
-0
src/metric.h less more
3535
3636 /*< public >*/
3737 GfsVariable * h[FTT_NEIGHBORS], * a;
38 gint level;
3839 };
3940
4041 #define GFS_METRIC_CUBED(obj) GTS_OBJECT_CAST (obj,\
+3
-0
src/variable.c less more
662662
663663 static void variable_stream_function_coarse_fine (FttCell * parent, GfsVariable * v)
664664 {
665 if (GFS_CELL_IS_BOUNDARY (parent))
666 return;
667
665668 GfsFunction * f = GFS_VARIABLE_FUNCTION (v)->f;
666669 FttCellChildren child;
667670 ftt_cell_children (parent, &child);