[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:25:17 UTC 2009
The following commit has been merged in the upstream branch:
commit cc79ac8a1baa9f8663af1931cdec3d3dce468dd9
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Oct 14 14:31:18 2009 +1100
"Cubed sphere" metric
darcs-hash:20091014033118-d4795-19d2084c7a7bfcdf59a4aae5dadaa3a1e18c738d.gz
diff --git a/src/Makefile.am b/src/Makefile.am
index 6b0aa37..ed8389e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,6 +59,7 @@ GFS_HDS = \
river.h \
moving.h \
balance.h \
+ cubed.h \
version.h
pkginclude_HEADERS = \
@@ -104,6 +105,7 @@ SRC = \
river.c \
moving.c \
balance.c \
+ cubed.c \
$(GFS_HDS)
domain.c: version.h
diff --git a/src/cubed.c b/src/cubed.c
new file mode 100644
index 0000000..17dd764
--- /dev/null
+++ b/src/cubed.c
@@ -0,0 +1,390 @@
+/* Gerris - The GNU Flow Solver
+ * Copyright (C) 2009 National Institute of Water and Atmospheric Research
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License as
+ * published by the Free Software Foundation; either version 2 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
+ * 02111-1307, USA.
+ */
+
+#include "cubed.h"
+#include <complex.h>
+#include "map.h"
+
+#define N 30
+
+#if 0
+/* Conformal mapping Taylor coefficients: from Rancic et al, 1996, Table B.1 */
+
+static double A[N] = {
+ 1.47713062600964, -0.38183510510174, -0.05573058001191, -0.01895883606818, -0.00791315785221,
+ -0.00486625437708, -0.00329251751279, -0.00235481488325, -0.00175870527475, -0.00135681133278,
+ -0.00107459847699, -0.00086944475948, -0.00071607115121, -0.00059867100093, -0.00050699063239,
+ -0.00043415191279, -0.00037541003286, -0.00032741060100, -0.00028773091482, -0.00025458777519,
+ -0.00022664642371, -0.00020289261022, -0.00018254510830, -0.00016499474461, -0.00014976117168,
+ -0.00013646173946, -0.00012478875823, -0.00011449267279, -0.00010536946150, -0.00009725109376
+};
+
+static double B[N] = {
+ 0.67698819751739, 0.11847293456554, 0.05317178134668, 0.02965810434052, 0.01912447304028,
+ 0.01342565621117, 0.00998873323180, 0.00774868996406, 0.00620346979888, 0.00509010874883,
+ 0.00425981184328, 0.00362308956077, 0.00312341468940, 0.00272360948942, 0.00239838086555,
+ 0.00213001905118, 0.00190581316131, 0.00171644156404, 0.00155493768255, 0.00141600715207,
+ 0.00129556597754, 0.00119042140226, 0.00109804711790, 0.00101642216628, 0.00094391366522,
+ 0.00087919021224, 0.00082115710311, 0.00076890728775, 0.00072168382969, 0.00067885087750
+};
+
+#else
+/* Conformal mapping Taylor coefficients: from map_xy2xyz.m from mitgcm */
+
+static double A[N] = {
+ 1.47713057321600, -0.38183513110512, -0.05573055466344, -0.01895884801823, -0.00791314396396,
+ -0.00486626515498, -0.00329250387158, -0.00235482619663, -0.00175869000970, -0.00135682443774,
+ -0.00107458043205, -0.00086946107050, -0.00071604933286, -0.00059869243613, -0.00050696402446,
+ -0.00043418115349, -0.00037537743098, -0.00032745130951, -0.00028769063795, -0.00025464473946,
+ -0.00022659577923, -0.00020297175587, -0.00018247947703, -0.00016510295548, -0.00014967258633,
+ -0.00013660647356, -0.00012466390509, -0.00011468147908, -0.00010518717478, -0.00009749136078
+};
+
+static double B[N] = {
+ 0.67698822171341, 0.11847295533659, 0.05317179075349, 0.02965811274764, 0.01912447871071,
+ 0.01342566129383, 0.00998873721022, 0.00774869352561, 0.00620347278164, 0.00509011141874,
+ 0.00425981415542, 0.00362309163280, 0.00312341651697, 0.00272361113245, 0.00239838233411,
+ 0.00213002038153, 0.00190581436893, 0.00171644267546, 0.00155493871562, 0.00141600812949,
+ 0.00129556691848, 0.00119042232809, 0.00109804804853, 0.00101642312253, 0.00094391466713,
+ 0.00087919127990, 0.00082115825576, 0.00076890854394, 0.00072168520663, 0.00067885239089
+};
+#endif
+
+static complex WofZ (complex Z)
+{
+ complex W = 0.;
+ int n = N;
+ while (n-- > 0)
+ W = (W + A[n])*Z;
+ return W;
+}
+
+static complex ZofW (complex W)
+{
+ complex Z = 0.;
+ int n = N;
+ while (n-- > 0)
+ Z = (Z + B[n])*W;
+ return Z;
+}
+
+/* I^(1/3) */
+#define I3 (0.86602540378444 + I/2.)
+/* sqrt (3.) - 1. */
+#define RA 0.73205080756888
+
+/* Conformal mapping of a cube onto a sphere. Maps (x,y) on the
+ * north-pole face of a cube to (X,Y,Z) coordinates in physical space.
+ *
+ * Based on f77 code from Jim Purser & Misha Rancic.
+ *
+ * Face is oriented normal to Z-axis with X and Y increasing with x
+ * and y.
+ *
+ * valid ranges: -1 < x < 1 -1 < y < 1.
+ */
+static void cmap_xy2XYZ (double x, double y, double * X, double * Y, double * Z)
+{
+ int kx = x < 0., ky = y < 0.;
+ x = fabs (x); y = fabs (y);
+ int kxy = y > x;
+
+ if (kxy) {
+ double tmp = x;
+ x = 1. - y;
+ y = 1. - tmp;
+ }
+ else {
+ x = 1. - x;
+ y = 1. - y;
+ }
+
+ complex z = (x + I*y)/2.;
+ complex W;
+ if (cabs (z) > 0.) {
+ z = z*z*z*z;
+ W = WofZ (z);
+ W = I3*cpow (W*I, 1./3.);
+ }
+ else
+ W = 0.;
+ complex cb = I - 1.;
+ complex cc = RA*cb/2.;
+ W = (W - RA)/(cb + cc*W);
+ *X = creal (W);
+ *Y = cimag (W);
+ double H = 2./(1. + (*X)*(*X) + (*Y)*(*Y));
+ *X *= H;
+ *Y *= H;
+ *Z = H - 1.;
+
+ if (kxy) {
+ double tmp = *X;
+ *X = *Y;
+ *Y = tmp;
+ }
+ if (kx)
+ *X = - *X;
+ if (ky)
+ *Y = - *Y;
+}
+
+/* Conformal mapping of a sphere onto a cube. Maps (X,Y,Z) coordinates
+ * in physical space to (x,y) on the north-pole face of a cube.
+ *
+ * This is the inverse transform of cmap_xy2XYZ().
+ */
+static void cmap_XYZ2xy (double X, double Y, double Z, double * x, double * y)
+{
+ int kx = X < 0., ky = Y < 0.;
+ X = fabs (X); Y = fabs (Y);
+ int kxy = Y > X;
+
+ if (kxy) {
+ double tmp = X;
+ X = Y;
+ Y = tmp;
+ }
+
+ double H = Z + 1.;
+ X /= H; Y /= H;
+ complex W = X + Y*I;
+ complex cb = I - 1.;
+ complex cc = RA*cb/2.;
+ W = (W*cb + RA)/(1. - W*cc);
+ W = W/I3;
+ W = W*W*W;
+ W /= I;
+ complex z = ZofW (W);
+ z = cpow (z, 1./4.)*2.;
+ *x = fabs (creal (z));
+ *y = fabs (cimag (z));
+
+ if (kxy) {
+ *x = 1. - *x;
+ *y = 1. - *y;
+ }
+ else {
+ double tmp = *x;
+ *x = 1. - *y;
+ *y = 1. - tmp;
+ }
+ if (kx)
+ *x = - *x;
+ if (ky)
+ *y = - *y;
+}
+
+static double angle_between_vectors (const double v1[3], const double v2[3])
+{
+ return acos (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
+}
+
+static void plane_normal (const double p1[3], const double p2[3], double plane[3])
+{
+ plane[0] = p1[1]*p2[2] - p1[2]*p2[1];
+ plane[1] = p1[2]*p2[0] - p1[0]*p2[2];
+ plane[2] = p1[0]*p2[1] - p1[1]*p2[0];
+ double mag = sqrt (plane[0]*plane[0] + plane[1]*plane[1] + plane[2]*plane[2]);
+ plane[0] /= mag;
+ plane[1] /= mag;
+ plane[2] /= mag;
+}
+
+static double excess_of_quad (const double v1[3], const double v2[3],
+ const double v3[3], const double v4[3])
+{
+ double plane1[3], plane2[3], plane3[3], plane4[3];
+
+ plane_normal (v1, v2, plane1);
+ plane_normal (v2, v3, plane2);
+ plane_normal (v3, v4, plane3);
+ plane_normal (v4, v1, plane4);
+
+ return 2.*M_PI -
+ angle_between_vectors (plane2, plane1) -
+ angle_between_vectors (plane3, plane2) -
+ angle_between_vectors (plane4, plane3) -
+ angle_between_vectors (plane1, plane4);
+}
+
+/* GfsMapCubed: Object */
+
+static void gfs_map_cubed_read (GtsObject ** o, GtsFile * fp)
+{
+ /* this mapping cannot be used independently from GfsMetricCubed */
+}
+
+static void gfs_map_cubed_write (GtsObject * o, FILE * fp)
+{
+ /* this mapping cannot be used independently from GfsMetricCubed */
+}
+
+static void gfs_map_cubed_class_init (GfsMapClass * klass)
+{
+ GTS_OBJECT_CLASS (klass)->read = gfs_map_cubed_read;
+ GTS_OBJECT_CLASS (klass)->write = gfs_map_cubed_write;
+}
+
+static void map_cubed_transform (GfsMap * map, const FttVector * src, FttVector * dest)
+{
+ GfsSimulation * sim = gfs_object_simulation (map);
+ double lon = src->x*M_PI/180., lat = src->y*M_PI/180.;
+ double X = cos (lat)*sin (lon), Y = sin (lat), Z = sqrt (1. - X*X - Y*Y);
+ double x, y;
+ cmap_XYZ2xy (X, Y, Z, &x, &y);
+ dest->x = x/2.*sim->physical_params.L;
+ dest->y = y/2.*sim->physical_params.L;
+ dest->z = src->z;
+}
+
+static void map_cubed_inverse (GfsMap * map, const FttVector * src, FttVector * dest)
+{
+ GfsSimulation * sim = gfs_object_simulation (map);
+ double X, Y, Z;
+ cmap_xy2XYZ (2.*src->x/sim->physical_params.L, 2.*src->y/sim->physical_params.L, &X, &Y, &Z);
+ double lat = asin (Y);
+ double lon = asin (X/cos (lat));
+ dest->x = lon*180./M_PI;
+ dest->y = lat*180./M_PI;
+ dest->z = src->z;
+}
+
+static void gfs_map_cubed_init (GfsMap * map)
+{
+ map->transform = map_cubed_transform;
+ map->inverse = map_cubed_inverse;
+}
+
+static GfsMapClass * gfs_map_cubed_class (void)
+{
+ static GfsMapClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_map_cubed_info = {
+ "GfsMapCubed",
+ sizeof (GfsMap),
+ sizeof (GfsMapClass),
+ (GtsObjectClassInitFunc) gfs_map_cubed_class_init,
+ (GtsObjectInitFunc) gfs_map_cubed_init,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_map_class ()), &gfs_map_cubed_info);
+ }
+
+ return klass;
+}
+
+/* GfsMetricCubed: Object */
+
+static gdouble metric_cubed_face_metric (const GfsDomain * domain, const FttCellFace * face,
+ const gpointer data)
+{
+ if (face->d/2 > 1)
+ return 1.;
+ FttVector p;
+ ftt_face_pos (face, &p);
+ double h = ftt_cell_size (face->cell);
+ double v1[3], v2[3];
+ v1[0] = v2[0] = 2.*p.x; v1[1] = v2[1] = 2.*p.y;
+ FttComponent c = FTT_ORTHOGONAL_COMPONENT (face->d/2);
+ v1[c] -= h; v2[c] += h;
+ cmap_xy2XYZ (v1[0], v1[1], &v1[0], &v1[1], &v1[2]);
+ cmap_xy2XYZ (v2[0], v2[1], &v2[0], &v2[1], &v2[2]);
+ return 2.*angle_between_vectors (v1, v2)/(M_PI*h);
+}
+
+static gdouble metric_cubed_cell_metric (const GfsDomain * domain, const FttCell * cell,
+ const gpointer data)
+{
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ double h = ftt_cell_size (cell);
+ double v1[3], v2[3], v3[3], v4[3];
+ v1[0] = 2.*p.x - h; v1[1] = 2.*p.y - h;
+ cmap_xy2XYZ (v1[0], v1[1], &v1[0], &v1[1], &v1[2]);
+ v2[0] = 2.*p.x - h; v2[1] = 2.*p.y + h;
+ cmap_xy2XYZ (v2[0], v2[1], &v2[0], &v2[1], &v2[2]);
+ v3[0] = 2.*p.x + h; v3[1] = 2.*p.y + h;
+ cmap_xy2XYZ (v3[0], v3[1], &v3[0], &v3[1], &v3[2]);
+ v4[0] = 2.*p.x + h; v4[1] = 2.*p.y - h;
+ cmap_xy2XYZ (v4[0], v4[1], &v4[0], &v4[1], &v4[2]);
+ return 4.*excess_of_quad (v1, v2, v3, v4)/(M_PI*M_PI*h*h);
+}
+
+static gdouble metric_cubed_solid_metric (const GfsDomain * domain, const FttCell * cell,
+ const gpointer data)
+{
+ g_assert (GFS_IS_MIXED (cell));
+ g_assert_not_implemented ();
+ return 1.;
+}
+
+static void metric_cubed_read (GtsObject ** o, GtsFile * fp)
+{
+ (* GTS_OBJECT_CLASS (gfs_metric_cubed_class ())->parent_class->read) (o, fp);
+ if (fp->type == GTS_ERROR)
+ return;
+
+ GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
+ if (domain->metric_data || domain->face_metric || domain->cell_metric || domain->solid_metric) {
+ gts_file_error (fp, "cannot use multiple metrics (yet)");
+ return;
+ }
+
+ GtsObject * map = gts_object_new (GTS_OBJECT_CLASS (gfs_map_cubed_class ()));
+ gfs_object_simulation_set (map, domain);
+ gts_container_add (GTS_CONTAINER (GFS_SIMULATION (domain)->maps), GTS_CONTAINEE (map));
+
+ domain->face_metric = metric_cubed_face_metric;
+ domain->cell_metric = metric_cubed_cell_metric;
+ domain->solid_metric = metric_cubed_solid_metric;
+}
+
+static void metric_cubed_class_init (GtsObjectClass * klass)
+{
+ klass->read = metric_cubed_read;
+}
+
+static void metric_cubed_init (GfsEvent * m)
+{
+ m->istep = G_MAXINT/2;
+}
+
+GfsEventClass * gfs_metric_cubed_class (void)
+{
+ static GfsEventClass * klass = NULL;
+
+ if (klass == NULL) {
+ GtsObjectClassInfo gfs_metric_cubed_info = {
+ "GfsMetricCubed",
+ sizeof (GfsEvent),
+ sizeof (GfsEventClass),
+ (GtsObjectClassInitFunc) metric_cubed_class_init,
+ (GtsObjectInitFunc) metric_cubed_init,
+ (GtsArgSetFunc) NULL,
+ (GtsArgGetFunc) NULL
+ };
+ klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_event_class ()),
+ &gfs_metric_cubed_info);
+ }
+
+ return klass;
+}
diff --git a/src/init.h b/src/cubed.h
similarity index 72%
copy from src/init.h
copy to src/cubed.h
index 425b2a2..f0e413c 100644
--- a/src/init.h
+++ b/src/cubed.h
@@ -1,5 +1,5 @@
/* Gerris - The GNU Flow Solver
- * Copyright (C) 2001 National Institute of Water and Atmospheric Research
+ * Copyright (C) 2009 National Institute of Water and Atmospheric Research
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
@@ -17,21 +17,24 @@
* 02111-1307, USA.
*/
-#ifndef __INIT_H__
-#define __INIT_H__
-
-#include <gts.h>
+#ifndef __CUBED_H__
+#define __CUBED_H__
#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */
-GtsObjectClass ** gfs_classes (void);
-void gfs_init (int * argc,
- char *** argv);
+#include "event.h"
+
+/* GfsMetricCubed: Header */
+
+#define GFS_IS_METRIC_CUBED(obj) (gts_object_is_from_class (obj,\
+ gfs_metric_cubed_class ()))
+
+GfsEventClass * gfs_metric_cubed_class (void);
#ifdef __cplusplus
}
#endif /* __cplusplus */
-#endif /* __INIT_H__ */
+#endif /* __CUBED_H__ */
diff --git a/src/init.c b/src/init.c
index d42ccdf..425957e 100644
--- a/src/init.c
+++ b/src/init.c
@@ -44,6 +44,7 @@
#include "river.h"
#include "balance.h"
#include "map.h"
+#include "cubed.h"
#include "modules.h"
@@ -166,6 +167,7 @@ GtsObjectClass ** gfs_classes (void)
gfs_init_wave_class (),
gfs_metric_lon_lat_class (),
+ gfs_metric_cubed_class (),
gfs_adapt_class (),
gfs_adapt_vorticity_class (),
diff --git a/src/surface.h b/src/surface.h
index 430ae11..d6fdae1 100644
--- a/src/surface.h
+++ b/src/surface.h
@@ -20,6 +20,10 @@
#ifndef __SURFACE_H__
#define __SURFACE_H__
+#ifdef __cplusplus
+extern "C" {
+#endif /* __cplusplus */
+
#include <gts.h>
#include "ftt.h"
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list