[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8

Stephane Popinet popinet at users.sf.net
Tue Nov 24 12:25:18 UTC 2009


The following commit has been merged in the upstream branch:
commit 4623e3bd3f4effbaf94af3cbf18de2743ab5e122
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Fri Oct 16 14:57:51 2009 +1100

    Optimisation of GfsMetricLonLat and GfsMetricCubed
    
    darcs-hash:20091016035751-d4795-ac1c93001f065e07aa7939926070492fc3053ef9.gz

diff --git a/src/Makefile.am b/src/Makefile.am
index ed8389e..d0ec605 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,7 +59,7 @@ GFS_HDS = \
 	river.h \
 	moving.h \
 	balance.h \
-	cubed.h \
+	metric.h \
 	version.h
 
 pkginclude_HEADERS = \
@@ -105,7 +105,7 @@ SRC = \
 	river.c \
 	moving.c \
 	balance.c \
-	cubed.c \
+	metric.c \
 	$(GFS_HDS)
 
 domain.c: version.h
diff --git a/src/cubed.c b/src/cubed.c
deleted file mode 100644
index 17dd764..0000000
--- a/src/cubed.c
+++ /dev/null
@@ -1,390 +0,0 @@
-/* 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/cubed.h b/src/cubed.h
deleted file mode 100644
index f0e413c..0000000
--- a/src/cubed.h
+++ /dev/null
@@ -1,40 +0,0 @@
-/* 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.  
- */
-
-#ifndef __CUBED_H__
-#define __CUBED_H__
-
-#ifdef __cplusplus
-extern "C" {
-#endif /* __cplusplus */
-
-#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 /* __CUBED_H__ */
diff --git a/src/init.c b/src/init.c
index 425957e..38d4a4e 100644
--- a/src/init.c
+++ b/src/init.c
@@ -44,7 +44,7 @@
 #include "river.h"
 #include "balance.h"
 #include "map.h"
-#include "cubed.h"
+#include "metric.h"
 
 #include "modules.h"
 
diff --git a/src/metric.c b/src/metric.c
new file mode 100644
index 0000000..40c18b0
--- /dev/null
+++ b/src/metric.c
@@ -0,0 +1,693 @@
+/* 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 "metric.h"
+#include <complex.h>
+#include "map.h"
+#include "solid.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.;
+  return GFS_VALUE (face->cell, GFS_METRIC_CUBED (data)->h[face->d]);
+}
+
+static gdouble metric_cubed_cell_metric (const GfsDomain * domain, const FttCell * cell, 
+					 const gpointer data)
+{
+  return GFS_VALUE (cell, GFS_METRIC_CUBED (data)->a);
+}
+
+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 none (FttCell * parent, GfsVariable * v)
+{
+}
+
+static void cubed_coarse_fine (FttCell * parent, GfsVariable * a)
+{
+  GfsMetricCubed * cubed = GTS_OBJECT (a)->reserved;
+  FttCellChildren child;
+  ftt_cell_children (parent, &child);
+  FttVector p;
+  ftt_cell_pos (parent, &p);
+  p.x *= 2.; p.y *= 2.;
+  double h = ftt_cell_size (parent);
+  double v0[3];
+  v0[0] = p.x;     v0[1] = p.y;
+  cmap_xy2XYZ (v0[0], v0[1], &v0[0], &v0[1], &v0[2]);
+  double v1[3], v2[3], v3[3], v4[3];
+  v1[0] = p.x - h; v1[1] = p.y - h;
+  cmap_xy2XYZ (v1[0], v1[1], &v1[0], &v1[1], &v1[2]);
+  v2[0] = p.x - h; v2[1] = p.y + h;
+  cmap_xy2XYZ (v2[0], v2[1], &v2[0], &v2[1], &v2[2]);
+  v3[0] = p.x + h; v3[1] = p.y + h;
+  cmap_xy2XYZ (v3[0], v3[1], &v3[0], &v3[1], &v3[2]);
+  v4[0] = p.x + h; v4[1] = p.y - h;
+  cmap_xy2XYZ (v4[0], v4[1], &v4[0], &v4[1], &v4[2]);
+  double v5[3], v6[3], v7[3], v8[3];
+  v5[0] = p.x; v5[1] = p.y - h;
+  cmap_xy2XYZ (v5[0], v5[1], &v5[0], &v5[1], &v5[2]);
+  v6[0] = p.x - h; v6[1] = p.y;
+  cmap_xy2XYZ (v6[0], v6[1], &v6[0], &v6[1], &v6[2]);
+  v7[0] = p.x; v7[1] = p.y + h;
+  cmap_xy2XYZ (v7[0], v7[1], &v7[0], &v7[1], &v7[2]);
+  v8[0] = p.x + h; v8[1] = p.y;
+  cmap_xy2XYZ (v8[0], v8[1], &v8[0], &v8[1], &v8[2]);
+
+  GFS_VALUE (child.c[0], a) = 16.*excess_of_quad (v6, v0, v7, v2)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[1], a) = 16.*excess_of_quad (v0, v8, v3, v7)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[2], a) = 16.*excess_of_quad (v1, v5, v0, v6)/(M_PI*M_PI*h*h);
+  GFS_VALUE (child.c[3], a) = 16.*excess_of_quad (v5, v4, v8, v0)/(M_PI*M_PI*h*h);
+
+  GFS_VALUE (child.c[0], cubed->h[0]) = GFS_VALUE (child.c[1], cubed->h[1]) = 
+    4.*angle_between_vectors (v0, v7)/(M_PI*h);
+  GFS_VALUE (child.c[0], cubed->h[3]) = GFS_VALUE (child.c[2], cubed->h[2]) = 
+    4.*angle_between_vectors (v0, v6)/(M_PI*h);
+  GFS_VALUE (child.c[2], cubed->h[0]) = GFS_VALUE (child.c[3], cubed->h[1]) = 
+    4.*angle_between_vectors (v0, v5)/(M_PI*h);
+  GFS_VALUE (child.c[1], cubed->h[3]) = GFS_VALUE (child.c[3], cubed->h[2]) = 
+    4.*angle_between_vectors (v0, v8)/(M_PI*h);
+
+  GFS_VALUE (child.c[0], cubed->h[2]) = 4.*angle_between_vectors (v2, v7)/(M_PI*h);
+  GFS_VALUE (child.c[0], cubed->h[1]) = 4.*angle_between_vectors (v2, v6)/(M_PI*h);
+  GFS_VALUE (child.c[1], cubed->h[2]) = 4.*angle_between_vectors (v3, v7)/(M_PI*h);
+  GFS_VALUE (child.c[1], cubed->h[0]) = 4.*angle_between_vectors (v3, v8)/(M_PI*h);
+  GFS_VALUE (child.c[2], cubed->h[3]) = 4.*angle_between_vectors (v1, v5)/(M_PI*h);
+  GFS_VALUE (child.c[2], cubed->h[1]) = 4.*angle_between_vectors (v1, v6)/(M_PI*h);
+  GFS_VALUE (child.c[3], cubed->h[0]) = 4.*angle_between_vectors (v4, v8)/(M_PI*h);
+  GFS_VALUE (child.c[3], cubed->h[3]) = 4.*angle_between_vectors (v4, v5)/(M_PI*h);
+}
+
+static void cubed_fine_coarse (FttCell * parent, GfsVariable * a)
+{
+  GfsMetricCubed * cubed = GTS_OBJECT (a)->reserved;
+  FttCellChildren child;
+  guint n;
+
+  ftt_cell_children (parent, &child);
+  gdouble va = 0.;
+  for (n = 0; n < 4; n++)
+    va += GFS_VALUE (child.c[n], a);
+  GFS_VALUE (parent, a) = va/4;
+
+  GFS_VALUE (parent, cubed->h[0]) = (GFS_VALUE (child.c[1], cubed->h[0]) +
+				     GFS_VALUE (child.c[3], cubed->h[0]))/2.;
+  GFS_VALUE (parent, cubed->h[1]) = (GFS_VALUE (child.c[0], cubed->h[1]) +
+				     GFS_VALUE (child.c[2], cubed->h[1]))/2.;
+  GFS_VALUE (parent, cubed->h[2]) = (GFS_VALUE (child.c[0], cubed->h[2]) +
+				     GFS_VALUE (child.c[1], cubed->h[2]))/2.;
+  GFS_VALUE (parent, cubed->h[3]) = (GFS_VALUE (child.c[2], cubed->h[3]) +
+				     GFS_VALUE (child.c[3], cubed->h[3]))/2.;
+}
+
+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;
+  }
+
+  GfsMetricCubed * cubed = GFS_METRIC_CUBED (*o);
+  FttDirection d;
+  for (d = 0; d < FTT_NEIGHBORS; d++) {
+    gchar * name = g_strdup_printf ("Ch%d", d);
+    cubed->h[d] = gfs_domain_get_or_add_variable (domain, name, "Cubed face metric");
+    cubed->h[d]->fine_coarse = cubed->h[d]->coarse_fine = none;
+    g_free (name);
+  }
+  cubed->a = gfs_domain_get_or_add_variable (domain, "Ca", "Cubed cell metric");
+  GTS_OBJECT (cubed->a)->reserved = cubed;
+  cubed->a->coarse_fine = cubed_coarse_fine;
+  cubed->a->fine_coarse = cubed_fine_coarse;
+
+  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->metric_data  = cubed;
+  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 (GfsMetricCubed),
+      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;
+}
+
+/* GfsMapLonLat: Header */
+
+typedef struct _GfsMapLonLat         GfsMapLonLat;
+
+struct _GfsMapLonLat {
+  /*< private >*/
+  GfsMap parent;
+
+  /*< public >*/
+  gdouble r;
+};
+
+#define GFS_MAP_LONLAT(obj)            GTS_OBJECT_CAST (obj,\
+					         GfsMapLonLat,\
+					         gfs_map_lonlat_class ())
+#define GFS_IS_MAP_LONLAT(obj)         (gts_object_is_from_class (obj,\
+						 gfs_map_lonlat_class ()))
+
+static GfsMapClass * gfs_map_lonlat_class      (void);
+
+/* GfsMapLonLat: Object */
+
+static void gfs_map_lonlat_read (GtsObject ** o, GtsFile * fp)
+{
+  /* this mapping cannot be used independently from GfsMetricLonLat */
+}
+
+static void gfs_map_lonlat_write (GtsObject * o, FILE * fp)
+{
+  /* this mapping cannot be used independently from GfsMetricLonLat */
+}
+
+static void gfs_map_lonlat_class_init (GfsMapClass * klass)
+{
+  GTS_OBJECT_CLASS (klass)->read = gfs_map_lonlat_read;
+  GTS_OBJECT_CLASS (klass)->write = gfs_map_lonlat_write;
+}
+
+static void map_lonlat_transform (GfsMap * map, const FttVector * src, FttVector * dest)
+{
+  dest->x = src->x*M_PI/180.*GFS_MAP_LONLAT (map)->r;
+  dest->y = src->y*M_PI/180.*GFS_MAP_LONLAT (map)->r;
+  dest->z = src->z;
+}
+
+static void map_lonlat_inverse (GfsMap * map, const FttVector * src, FttVector * dest)
+{
+  dest->x = src->x*180./(M_PI*GFS_MAP_LONLAT (map)->r);
+  dest->y = src->y*180./(M_PI*GFS_MAP_LONLAT (map)->r);
+  dest->z = src->z;
+}
+
+static void gfs_map_lonlat_init (GfsMap * map)
+{
+  map->transform = map_lonlat_transform;
+  map->inverse =   map_lonlat_inverse;
+  GFS_MAP_LONLAT (map)->r = 1.;
+}
+
+static GfsMapClass * gfs_map_lonlat_class (void)
+{
+  static GfsMapClass * klass = NULL;
+
+  if (klass == NULL) {
+    GtsObjectClassInfo gfs_map_lonlat_info = {
+      "GfsMapLonLat",
+      sizeof (GfsMapLonLat),
+      sizeof (GfsMapClass),
+      (GtsObjectClassInitFunc) gfs_map_lonlat_class_init,
+      (GtsObjectInitFunc) gfs_map_lonlat_init,
+      (GtsArgSetFunc) NULL,
+      (GtsArgGetFunc) NULL
+    };
+    klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_map_class ()), &gfs_map_lonlat_info);
+  }
+
+  return klass;
+}
+
+/* GfsMetricLonLat: Object */
+
+static void metric_lon_lat_write (GtsObject * o, FILE * fp)
+{
+  (* GTS_OBJECT_CLASS (gfs_metric_lon_lat_class ())->parent_class->write) (o, fp);
+  fprintf (fp, " %g", GFS_METRIC_LON_LAT (o)->r);
+}
+
+static gdouble metric_lon_lat_face_metric (const GfsDomain * domain, const FttCellFace * face, 
+					   const gpointer data)
+{
+  if (face->d/2 != FTT_Y)
+    return 1.;
+  return face->d == 2 ? 
+    GFS_VALUE (face->cell, GFS_METRIC_LON_LAT (data)->h2) :
+    GFS_VALUE (face->cell, GFS_METRIC_LON_LAT (data)->h3);
+}
+
+static gdouble metric_lon_lat_cell_metric (const GfsDomain * domain, const FttCell * cell, 
+					   const gpointer data)
+{
+  return GFS_VALUE (cell, GFS_METRIC_LON_LAT (data)->a);
+}
+
+static gdouble metric_lon_lat_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 lonlat_coarse_fine (FttCell * parent, GfsVariable * a)
+{
+  GfsMetricLonLat * lonlat = GTS_OBJECT (a)->reserved;
+  FttCellChildren child;
+  ftt_cell_children (parent, &child);
+  FttVector p;
+  ftt_cell_pos (parent, &p);
+  double theta0 = gfs_object_simulation (lonlat)->physical_params.L/lonlat->r;
+  double theta = p.y*theta0;
+  double h = ftt_cell_size (parent);
+  double dtheta = h*theta0/2.;
+  double theta1 = theta + dtheta;
+  double theta2 = theta - dtheta;
+  double sintheta = sin (theta);
+
+  GFS_VALUE (child.c[0], a) = GFS_VALUE (child.c[1], a) = (sin (theta1) - sintheta)/dtheta;
+  GFS_VALUE (child.c[2], a) = GFS_VALUE (child.c[3], a) = (sintheta - sin (theta2))/dtheta;
+
+  GFS_VALUE (child.c[0], lonlat->h2) = GFS_VALUE (child.c[1], lonlat->h2) = cos (theta1);
+  GFS_VALUE (child.c[0], lonlat->h3) = GFS_VALUE (child.c[1], lonlat->h3) = 
+    GFS_VALUE (child.c[2], lonlat->h2) = GFS_VALUE (child.c[3], lonlat->h2) = 
+    cos (theta);
+  GFS_VALUE (child.c[2], lonlat->h3) = GFS_VALUE (child.c[3], lonlat->h3) = cos (theta2);
+}
+
+static void lonlat_fine_coarse (FttCell * parent, GfsVariable * a)
+{
+  GfsMetricLonLat * lonlat = GTS_OBJECT (a)->reserved;
+  FttCellChildren child;
+  guint n;
+
+  ftt_cell_children (parent, &child);
+  gdouble va = 0.;
+  for (n = 0; n < 4; n++)
+    va += GFS_VALUE (child.c[n], a);
+  GFS_VALUE (parent, a) = va/4;
+
+  GFS_VALUE (parent, lonlat->h2) = (GFS_VALUE (child.c[0], lonlat->h2) +
+				    GFS_VALUE (child.c[1], lonlat->h2))/2.;
+  GFS_VALUE (parent, lonlat->h3) = (GFS_VALUE (child.c[3], lonlat->h3) +
+				    GFS_VALUE (child.c[2], lonlat->h3))/2.;
+}
+
+static void metric_lon_lat_read (GtsObject ** o, GtsFile * fp)
+{
+  (* GTS_OBJECT_CLASS (gfs_metric_lon_lat_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;
+  }
+
+  GFS_METRIC_LON_LAT (*o)->r = gfs_read_constant (fp, gfs_object_simulation (*o));
+  if (fp->type == GTS_ERROR)
+    return;
+  if (GFS_METRIC_LON_LAT (*o)->r <= 0.) {
+    gts_file_error (fp, "radius must be strictly positive");
+    return;
+  }
+
+  GfsMetricLonLat * lonlat = GFS_METRIC_LON_LAT (*o);
+  lonlat->h2 = gfs_domain_get_or_add_variable (domain, "Lh2", "LonLat face metric");
+  lonlat->h2->coarse_fine = lonlat->h2->fine_coarse = none;
+  lonlat->h3 = gfs_domain_get_or_add_variable (domain, "Lh3", "LonLat face metric");
+  lonlat->h3->coarse_fine = lonlat->h3->fine_coarse = none;
+  lonlat->a =  gfs_domain_get_or_add_variable (domain, "La",  "LonLat cell metric");
+  GTS_OBJECT (lonlat->a)->reserved = lonlat;
+  lonlat->a->coarse_fine = lonlat_coarse_fine;
+  lonlat->a->fine_coarse = lonlat_fine_coarse;
+
+  GtsObject * map = gts_object_new (GTS_OBJECT_CLASS (gfs_map_lonlat_class ()));
+  gfs_object_simulation_set (map, domain);
+  gts_container_add (GTS_CONTAINER (GFS_SIMULATION (domain)->maps), GTS_CONTAINEE (map));
+  GFS_MAP_LONLAT (map)->r = GFS_METRIC_LON_LAT (*o)->r;
+
+  domain->metric_data = *o;
+  domain->face_metric  = metric_lon_lat_face_metric;
+  domain->cell_metric  = metric_lon_lat_cell_metric;
+  domain->solid_metric = metric_lon_lat_solid_metric;
+}
+
+static void metric_lon_lat_class_init (GtsObjectClass * klass)
+{
+  klass->read = metric_lon_lat_read;
+  klass->write = metric_lon_lat_write;
+}
+
+static void metric_lon_lat_init (GfsMetricLonLat * m)
+{
+  GFS_EVENT (m)->istep = G_MAXINT/2;
+  m->r = 1.;
+}
+
+GfsEventClass * gfs_metric_lon_lat_class (void)
+{
+  static GfsEventClass * klass = NULL;
+
+  if (klass == NULL) {
+    GtsObjectClassInfo gfs_metric_lon_lat_info = {
+      "GfsMetricLonLat",
+      sizeof (GfsMetricLonLat),
+      sizeof (GfsEventClass),
+      (GtsObjectClassInitFunc) metric_lon_lat_class_init,
+      (GtsObjectInitFunc) metric_lon_lat_init,
+      (GtsArgSetFunc) NULL,
+      (GtsArgGetFunc) NULL
+    };
+    klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_event_class ()), 
+				  &gfs_metric_lon_lat_info);
+  }
+
+  return klass;
+}
diff --git a/src/metric.h b/src/metric.h
new file mode 100644
index 0000000..76843e1
--- /dev/null
+++ b/src/metric.h
@@ -0,0 +1,74 @@
+/* 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.  
+ */
+
+#ifndef __METRIC_H__
+#define __METRIC_H__
+
+#ifdef __cplusplus
+extern "C" {
+#endif /* __cplusplus */
+
+#include "event.h"
+
+/* GfsMetricCubed: Header */
+
+typedef struct _GfsMetricCubed GfsMetricCubed;
+
+struct _GfsMetricCubed {
+  /*< private >*/
+  GfsEvent parent;
+
+  /*< public >*/
+  GfsVariable * h[FTT_NEIGHBORS], * a;
+};
+
+#define GFS_METRIC_CUBED(obj)            GTS_OBJECT_CAST (obj,\
+					           GfsMetricCubed,\
+					           gfs_metric_cubed_class ())
+#define GFS_IS_METRIC_CUBED(obj)         (gts_object_is_from_class (obj,\
+						   gfs_metric_cubed_class ()))
+
+GfsEventClass * gfs_metric_cubed_class  (void);
+
+/* GfsMetricLonLat: Header */
+
+typedef struct _GfsMetricLonLat GfsMetricLonLat;
+
+struct _GfsMetricLonLat {
+  /*< private >*/
+  GfsEvent parent;
+
+  /*< public >*/
+  GfsVariable * a, * h2, * h3;
+  gdouble r;
+};
+
+#define GFS_METRIC_LON_LAT(obj)            GTS_OBJECT_CAST (obj,\
+					           GfsMetricLonLat,\
+					           gfs_metric_lon_lat_class ())
+#define GFS_IS_METRIC_LON_LAT(obj)         (gts_object_is_from_class (obj,\
+						   gfs_metric_lon_lat_class ()))
+
+GfsEventClass * gfs_metric_lon_lat_class  (void);
+
+#ifdef __cplusplus
+}
+#endif /* __cplusplus */
+
+#endif /* __METRIC_H__ */
diff --git a/src/simulation.c b/src/simulation.c
index 0b5b80c..34c2dac 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -2020,96 +2020,3 @@ GfsSimulationClass * gfs_axi_class (void)
 
   return klass;
 }
-
-/* GfsMetricLonLat: Object */
-
-static void metric_lon_lat_write (GtsObject * o, FILE * fp)
-{
-  (* GTS_OBJECT_CLASS (gfs_metric_lon_lat_class ())->parent_class->write) (o, fp);
-  fprintf (fp, " %g", GFS_METRIC_LON_LAT (o)->r);
-}
-
-static gdouble metric_lon_lat_face_metric (const GfsDomain * domain, const FttCellFace * face, 
-					   const gpointer data)
-{
-  if (face->d/2 != FTT_Y)
-    return 1.;
-  FttVector p;
-  ftt_face_pos (face, &p);
-  return cos (p.y*GFS_SIMULATION (domain)->physical_params.L/GFS_METRIC_LON_LAT (data)->r);
-}
-
-static gdouble metric_lon_lat_cell_metric (const GfsDomain * domain, const FttCell * cell, 
-					   const gpointer data)
-{
-  FttVector p;
-  gfs_cell_cm (cell, &p);
-  return cos (p.y*GFS_SIMULATION (domain)->physical_params.L/GFS_METRIC_LON_LAT (data)->r);
-}
-
-static gdouble metric_lon_lat_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_lon_lat_read (GtsObject ** o, GtsFile * fp)
-{
-  (* GTS_OBJECT_CLASS (gfs_metric_lon_lat_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;
-  }
-
-  GFS_METRIC_LON_LAT (*o)->r = gfs_read_constant (fp, gfs_object_simulation (*o));
-  if (fp->type == GTS_ERROR)
-    return;
-  if (GFS_METRIC_LON_LAT (*o)->r <= 0.) {
-    gts_file_error (fp, "radius must be strictly positive");
-    return;
-  }
-
-  domain->metric_data = *o;
-  domain->face_metric  = metric_lon_lat_face_metric;
-  domain->cell_metric  = metric_lon_lat_cell_metric;
-  domain->solid_metric = metric_lon_lat_solid_metric;
-}
-
-static void metric_lon_lat_class_init (GtsObjectClass * klass)
-{
-  klass->read = metric_lon_lat_read;
-  klass->write = metric_lon_lat_write;
-}
-
-static void metric_lon_lat_init (GfsMetricLonLat * m)
-{
-  GFS_EVENT (m)->istep = G_MAXINT/2;
-  m->r = 1.;
-}
-
-GfsEventClass * gfs_metric_lon_lat_class (void)
-{
-  static GfsEventClass * klass = NULL;
-
-  if (klass == NULL) {
-    GtsObjectClassInfo gfs_metric_lon_lat_info = {
-      "GfsMetricLonLat",
-      sizeof (GfsMetricLonLat),
-      sizeof (GfsEventClass),
-      (GtsObjectClassInitFunc) metric_lon_lat_class_init,
-      (GtsObjectInitFunc) metric_lon_lat_init,
-      (GtsArgSetFunc) NULL,
-      (GtsArgGetFunc) NULL
-    };
-    klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_event_class ()), 
-				  &gfs_metric_lon_lat_info);
-  }
-
-  return klass;
-}
diff --git a/src/simulation.h b/src/simulation.h
index 544ac95..f01c08d 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -147,26 +147,6 @@ GfsSimulationClass * gfs_poisson_class            (void);
 
 GfsSimulationClass * gfs_axi_class                (void);
 
-/* GfsMetricLonLat: Header */
-
-typedef struct _GfsMetricLonLat GfsMetricLonLat;
-
-struct _GfsMetricLonLat {
-  /*< private >*/
-  GfsEvent parent;
-
-  /*< public >*/
-  gdouble r;
-};
-
-#define GFS_METRIC_LON_LAT(obj)            GTS_OBJECT_CAST (obj,\
-					           GfsMetricLonLat,\
-					           gfs_metric_lon_lat_class ())
-#define GFS_IS_METRIC_LON_LAT(obj)         (gts_object_is_from_class (obj,\
-						   gfs_metric_lon_lat_class ()))
-
-GfsEventClass * gfs_metric_lon_lat_class  (void);
-
 #ifdef __cplusplus
 }
 #endif /* __cplusplus */
diff --git a/test/lonlat/coriolis/coriolis.gfs b/test/lonlat/coriolis/coriolis.gfs
index aa59a24..e334e54 100644
--- a/test/lonlat/coriolis/coriolis.gfs
+++ b/test/lonlat/coriolis/coriolis.gfs
@@ -33,10 +33,10 @@ Define LENGTH (150./180.*M_PI)
     PhysicalParams { L = LENGTH }
     MetricLonLat 1.
     Refine 8
-    InitFraction P (0.2 - acos(cos(x)*cos(y)))
+    InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
     Init {} { P = 0.2 + 1.8*P/LENGTH }
     Time { end = 1.4 }
-    SourceCoriolis 10.*sin(y)
+    SourceCoriolis 10.*sin(y*M_PI/180.)
     OutputTime { istep = 10 } stderr
     OutputSimulation { step = 0.4 } sim-%g.gfs
     OutputSimulation { step = 0.4 } sim-%g.txt { variables = U,V,P format = text }
diff --git a/test/lonlat/lonlat.gfs b/test/lonlat/lonlat.gfs
index 0135b77..6ca9f85 100644
--- a/test/lonlat/lonlat.gfs
+++ b/test/lonlat/lonlat.gfs
@@ -53,7 +53,7 @@ Define LENGTH (150./180.*M_PI)
     PhysicalParams { L = LENGTH }
     MetricLonLat 1.
     Refine 8
-    InitFraction P (0.2 - acos(cos(x)*cos(y)))
+    InitFraction P (0.2 - acos(cos(x*M_PI/180.)*cos(y*M_PI/180.)))
     Init {} { P = 0.2 + 1.8*P/LENGTH }
     Time { end = 0.9 }
     OutputTime { istep = 10 } stderr
@@ -62,9 +62,9 @@ Define LENGTH (150./180.*M_PI)
 #    OutputSimulation { istep = 10 } stdout
     EventScript { start = end } {
 	for i in 0.3 0.6 0.9; do
-	    if awk 'BEGIN { n1 = 0 }{
-              c = cos($1)*cos($2);
-              d = atan2(sqrt(1. - c*c),c)*180./3.14159265359
+	    if awk 'BEGIN { n1 = 0; pi = 3.14159265359 }{
+              c = cos($1*pi/180.)*cos($2*pi/180.);
+              d = atan2(sqrt(1. - c*c),c)*180./pi
               i = int(d*2.)
               x[i] += d
               y[i] += $6
@@ -92,8 +92,8 @@ Define LENGTH (150./180.*M_PI)
 	    fi
 
 	    cat <<EOF | gnuplot
-rdist(x,y)=acos(cos(x)*cos(y))*180./pi
-cdist(x,y)=sqrt(x*x+y*y)*180./pi
+rdist(x,y)=acos(cos(x*pi/180.)*cos(y*pi/180.))*180./pi
+cdist(x,y)=sqrt(x*x+y*y)
 set xlabel 'Angular distance (degree)'
 set ylabel 'Surface height'
 set xtics 0,22.5,90

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list