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

Stephane Popinet popinet at users.sf.net
Fri May 15 02:54:15 UTC 2009


The following commit has been merged in the upstream branch:
commit 931fcb58a93768489747a801434c978ebf40c39d
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Tue Apr 17 10:17:48 2007 +1000

    All surfaces are now defined using a new GfsSurface object
    
    Note that this is a new object, the previous "GfsSurface" has been replaced with
    "GfsSolid".
    
    darcs-hash:20070417001748-d4795-9faa02c177c31cbf6a1a040498021c5258d7e56e.gz

diff --git a/configure.in b/configure.in
index 1aca81c..b05959a 100644
--- a/configure.in
+++ b/configure.in
@@ -12,7 +12,7 @@ dnl are available for $ac_help expansion (don't we all *love* autoconf?)
 # set GFS_BINARY_AGE and GFS_INTERFACE_AGE to 0.
 #
 GFS_MAJOR_VERSION=1
-GFS_MINOR_VERSION=0
+GFS_MINOR_VERSION=1
 GFS_MICRO_VERSION=0
 GFS_INTERFACE_AGE=0
 GFS_BINARY_AGE=0
diff --git a/doc/examples/boussinesq/boussinesq.gfs b/doc/examples/boussinesq/boussinesq.gfs
index 401f9fe..cbb7321 100644
--- a/doc/examples/boussinesq/boussinesq.gfs
+++ b/doc/examples/boussinesq/boussinesq.gfs
@@ -36,7 +36,7 @@
   RefineSolid 8
 
   # Insert the solid boundary defined by cylinder.gts
-  Surface cylinder.gts
+  Solid cylinder.gts
 
   # Add a passive tracer called T
   VariableTracer {} T
diff --git a/doc/examples/cylinder/cylinder.gfs b/doc/examples/cylinder/cylinder.gfs
index e4905b0..3ae39c2 100644
--- a/doc/examples/cylinder/cylinder.gfs
+++ b/doc/examples/cylinder/cylinder.gfs
@@ -48,7 +48,7 @@
   Refine 6
 
   # Insert the solid boundary defined by cylinder.gts
-  Surface cylinder.gts
+  Solid cylinder.gts
 
   # Add a passive tracer called T
   VariableTracer {} T
diff --git a/doc/examples/cylinder/heated/heated.gfs b/doc/examples/cylinder/heated/heated.gfs
index 6d0a20e..8b78b31 100644
--- a/doc/examples/cylinder/heated/heated.gfs
+++ b/doc/examples/cylinder/heated/heated.gfs
@@ -34,7 +34,7 @@
   Refine 6
 
   # Insert the solid boundary defined by cylinder.gts
-  Surface cylinder.gts
+  Solid cylinder.gts
 
   # Add a passive tracer called T
   VariableTracer {} T
diff --git a/doc/examples/logo/logo.gfs b/doc/examples/logo/logo.gfs
index 7c25de5..b79bea2 100644
--- a/doc/examples/logo/logo.gfs
+++ b/doc/examples/logo/logo.gfs
@@ -42,7 +42,7 @@
     AdaptVorticity { istep = 1 } { cmax = 1e-2 maxlevel = 12 }
     OutputTime { istep = 1 } stderr
     OutputProjectionStats { istep = 1 } stderr
-    OutputSimulation { istep = 10 } stdout {}
+    OutputSimulation { istep = 10 } stdout
     OutputPPM { istep = 2 } { ppm2mpeg > logo.mpg } {
         v = Vorticity
         min = -0.1348 max = 6.22219
diff --git a/doc/examples/rt/rt.gfs b/doc/examples/rt/rt.gfs
index 9c78f68..12ce6b9 100644
--- a/doc/examples/rt/rt.gfs
+++ b/doc/examples/rt/rt.gfs
@@ -78,7 +78,7 @@
     min = 0 max = 1 v = T
   }
   OutputTiming { start = end } stderr
-  OutputSimulation { step = 0.1 } stdout {}
+  OutputSimulation { step = 0.1 } stdout
   EventScript { start = 0 } { echo "Save t-0.eps { format = EPS }" }
   EventScript { start = 0.7 step = 0.1 } { echo "Save t-$GfsTime.eps { format = EPS }" }
 }
diff --git a/doc/examples/tangaroa/tangaroa.gfs b/doc/examples/tangaroa/tangaroa.gfs
index 23fb2f7..3b8f76c 100644
--- a/doc/examples/tangaroa/tangaroa.gfs
+++ b/doc/examples/tangaroa/tangaroa.gfs
@@ -22,7 +22,7 @@
 #
 2 1 GfsSimulation GfsBox GfsGEdge {} {
   Time { end = 2 }
-  Surface tangaroa.gts
+  Solid tangaroa.gts
   Refine 5
   RefineSolid 9
   Init {} { U = 1. }
@@ -50,7 +50,7 @@
   EventSum { start = 1 istep = 1 } W*W SW2
 
   # Output simulation on standard output (to be read and displayed by GfsView)
-  OutputSimulation { istep = 4 } stdout {}
+  OutputSimulation { istep = 4 } stdout
   # Sends a command to GfsView to save a 1024x768 PPM image on standard output
   EventScript { istep = 4 } { echo "Save stdout { width = 1024 height = 768 }" }
 
diff --git a/doc/tutorial/tutorial.tex b/doc/tutorial/tutorial.tex
index 719711e..e5267b4 100644
--- a/doc/tutorial/tutorial.tex
+++ b/doc/tutorial/tutorial.tex
@@ -897,7 +897,7 @@ We can now insert this object in the simulation domain like this:
 \begin{verbatim}
 4 3 GfsSimulation GfsBox GfsGEdge {} {
   GfsTime { end = 0 }
-  GfsSurface half-cylinder.gts
+  GfsSolid half-cylinder.gts
 }
 GfsBox { left = GfsBoundaryInflowConstant 1 }
 GfsBox {}
@@ -912,7 +912,7 @@ add what mesh refinement we want and a few things to output:
 4 3 GfsSimulation GfsBox GfsGEdge {} {
   GfsTime { end = 9 }
   GfsRefine 6
-  GfsSurface half-cylinder.gts
+  GfsSolid half-cylinder.gts
   GfsInit {} { U = 1 }
   GfsOutputBoundaries {} boundaries
   GfsOutputTime { step = 0.02 } stdout
@@ -1009,10 +1009,9 @@ simulation from this point onward. If you edit it, you will see that
 the general structure is the same as usual but for five pretty big
 chunks of data. 
 
-The first chunk starts with {\tt GtsSurface} and is
+The first chunk starts with {\tt GfsSolid} and is
 just the data contained in {\tt half-cylinder.gts} but this time
-embedded directly (by using {\tt GtsSurface} rather than {\tt
-GfsSurface}) into the simulation file. The goal there is to have
+embedded directly into the simulation file. The goal there is to have
 fully self-contained simulations files which you can just move around
 without having to keep track of twenty different files.
 
@@ -1191,7 +1190,7 @@ Following this we can modify our simulation file:
 4 3 GfsSimulation GfsBox GfsGEdge {} {
   GfsTime { end = 9 }
   GfsRefine 7
-  GfsSurface half-cylinder.gts
+  GfsSolid half-cylinder.gts
   GfsInit {} { U = 1 }
 #  GfsOutputBoundaries {} boundaries
   GfsAdaptVorticity { istep = 1 } { maxlevel = 7 cmax = 1e-2 }
@@ -1343,7 +1342,7 @@ file like this:
 4 3 GfsSimulation GfsBox GfsGEdge {} {
   GfsTime { end = 9 }
   GfsRefine 7
-  GfsSurface half-cylinder.gts
+  GfsSolid half-cylinder.gts
   GfsVariableTracer {} T
   ...
   GfsOutputPPM { step = 0.02 } tracer.ppm {
diff --git a/src/Makefile.am b/src/Makefile.am
index 986ce28..cb4e371 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -54,6 +54,7 @@ GFS_HDS = \
 	levelset.h \
 	isocube.h \
 	cartesian.h \
+	surface.h \
 	version.h \
 	$(MPI_HDS)
 
@@ -91,6 +92,7 @@ SRC = \
 	myc.h \
 	myc2d.h \
 	cartesian.c \
+	surface.c \
 	$(GFS_HDS)
 
 simulation.c: version.h
diff --git a/src/domain.c b/src/domain.c
index 0628c74..2f1766c 100644
--- a/src/domain.c
+++ b/src/domain.c
@@ -917,21 +917,23 @@ void gfs_domain_traverse_mixed (GfsDomain * domain,
   gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_mixed, datum);
 }
 
-static void traverse_cut (GfsBox * box, gpointer * datum)
-{
-  FttCellTraverseCutFunc func = (FttCellTraverseCutFunc) datum[0];
-  gpointer data = datum[1];
-  FttTraverseType * order = datum[2];
-  FttTraverseFlags * flags = datum[3];
-  GtsSurface * s = datum[4];
+typedef struct {
+  FttCellTraverseCutFunc func;
+  gpointer data;
+  FttTraverseType order;
+  FttTraverseFlags flags;
+  GfsSurface * s;
+} TraverseCut;
 
-  gfs_cell_traverse_cut (box->root, s, *order, *flags, func, data);
+static void traverse_cut (GfsBox * box, TraverseCut * p)
+{
+  gfs_cell_traverse_cut (box->root, p->s, p->order, p->flags, p->func, p->data);
 }
 
 /**
  * gfs_domain_traverse_cut:
  * @domain: a #GfsDomain.
- * @s: a #GtsSurface.
+ * @s: a #GfsSurface.
  * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
  * %FTT_POST_ORDER. 
  * @flags: which types of children are to be visited.
@@ -941,42 +943,36 @@ static void traverse_cut (GfsBox * box, gpointer * datum)
  * Calls @func for each cell of @domain cut by @s.
  */
 void gfs_domain_traverse_cut (GfsDomain * domain,
-			      GtsSurface * s,
+			      GfsSurface * s,
 			      FttTraverseType order,
 			      FttTraverseFlags flags,
 			      FttCellTraverseCutFunc func,
 			      gpointer data)
 {
-  gpointer datum[5];
+  TraverseCut p;
 
-  datum[0] = func;
-  datum[1] = data;
-  datum[2] = &order;
-  datum[3] = &flags;
-  datum[4] = s;
+  p.func = func;
+  p.data= data;
+  p.order = order;
+  p.flags = flags;
+  p.s = s;
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (s != NULL);
   g_return_if_fail (func != NULL);
 
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_cut, datum);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_cut, &p);
 }
 
-static void traverse_cut_2D (GfsBox * box, gpointer * datum)
+static void traverse_cut_2D (GfsBox * box, TraverseCut * p)
 {
-  FttCellTraverseCutFunc func = (FttCellTraverseCutFunc) datum[0];
-  gpointer data = datum[1];
-  FttTraverseType * order = datum[2];
-  FttTraverseFlags * flags = datum[3];
-  GtsSurface * s = datum[4];
-
-  gfs_cell_traverse_cut_2D (box->root, s, *order, *flags, func, data);
+  gfs_cell_traverse_cut_2D (box->root, p->s, p->order, p->flags, p->func, p->data);
 }
 
 /**
  * gfs_domain_traverse_cut_2D:
  * @domain: a #GfsDomain.
- * @s: a #GtsSurface.
+ * @s: a #GfsSurface.
  * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
  * %FTT_POST_ORDER. 
  * @flags: which types of children are to be visited.
@@ -988,25 +984,25 @@ static void traverse_cut_2D (GfsBox * box, gpointer * datum)
  * The cells are flattened in the z-direction.
  */
 void gfs_domain_traverse_cut_2D (GfsDomain * domain,
-				 GtsSurface * s,
+				 GfsSurface * s,
 				 FttTraverseType order,
 				 FttTraverseFlags flags,
 				 FttCellTraverseCutFunc func,
 				 gpointer data)
 {
-  gpointer datum[5];
+  TraverseCut p;
 
-  datum[0] = func;
-  datum[1] = data;
-  datum[2] = &order;
-  datum[3] = &flags;
-  datum[4] = s;
+  p.func = func;
+  p.data = data;
+  p.order = order;
+  p.flags = flags;
+  p.s = s;
 
   g_return_if_fail (domain != NULL);
   g_return_if_fail (s != NULL);
   g_return_if_fail (func != NULL);
 
-  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_cut_2D, datum);
+  gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) traverse_cut_2D, &p);
 }
 
 static void box_depth (GfsBox * box, guint * depth)
diff --git a/src/domain.h b/src/domain.h
index a432d72..0055ee3 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -21,6 +21,7 @@
 #define __DOMAIN_H__
 
 #include "boundary.h"
+#include "surface.h"
 
 #ifdef __cplusplus
 extern "C" {
@@ -120,13 +121,13 @@ void         gfs_domain_traverse_mixed        (GfsDomain * domain,
 					       FttCellTraverseFunc func,
 					       gpointer data);
 void         gfs_domain_traverse_cut          (GfsDomain * domain,
-					       GtsSurface * s,
+					       GfsSurface * s,
 					       FttTraverseType order,
 					       FttTraverseFlags flags,
 					       FttCellTraverseCutFunc func,
 					       gpointer data);
 void         gfs_domain_traverse_cut_2D       (GfsDomain * domain,
-					       GtsSurface * s,
+					       GfsSurface * s,
 					       FttTraverseType order,
 					       FttTraverseFlags flags,
 					       FttCellTraverseCutFunc func,
diff --git a/src/event.c b/src/event.c
index cbfdc40..bc74d38 100644
--- a/src/event.c
+++ b/src/event.c
@@ -1616,8 +1616,7 @@ static void gfs_init_fraction_destroy (GtsObject * object)
 {
   GfsInitFraction * init = GFS_INIT_FRACTION (object);
 
-  if (init->surface)
-    gts_object_destroy (GTS_OBJECT (init->surface));
+  gts_object_destroy (GTS_OBJECT (init->surface));
 
   (* GTS_OBJECT_CLASS (gfs_init_fraction_class ())->parent_class->destroy) 
     (object);
@@ -1647,63 +1646,17 @@ static void gfs_init_fraction_read (GtsObject ** o, GtsFile * fp)
   }
   gts_file_next_token (fp);
 
-  if (fp->type != '{') {
-    FILE * f;
-    GtsFile * gf;
-
-    if (fp->type != GTS_STRING) {
-      gts_file_error (fp, "expecting a string (filename)\n");
-      return;
-    }
-    f = fopen (fp->token->str, "rt");
-    if (f == NULL) {
-      gts_file_error (fp, "cannot open file `%s'\n", fp->token->str);
-      return;
-    }
-    gf = gts_file_new (f);
-    if (gts_surface_read (init->surface, gf)) {
-      gts_file_error (fp, 
-		      "file `%s' is not a valid GTS file\n"
-		      "%s:%d:%d: %s",
-		      fp->token->str, fp->token->str,
-		      gf->line, gf->pos, gf->error);
-      gts_file_destroy (gf);
-      fclose (f);
-      return;
-    }
-    gts_file_destroy (gf);
-    fclose (f);
-  }
-  else { /* embedded GTS file */
-    fp->scope_max++;
-    gts_file_next_token (fp);
-    if (gts_surface_read (init->surface, fp))
-      return;
-    if (fp->type != '}') {
-      gts_file_error (fp, "expecting a closing brace");
-      return;
-    }
-    fp->scope_max--;
-  }
-  
-  if (!gts_surface_is_orientable (init->surface)) {
-    gts_file_error (fp, "surface is not orientable");
-    return;
-  }
-  
-  gts_file_next_token (fp);
+  gfs_surface_read (init->surface, fp);
 }
 
 static void gfs_init_fraction_write (GtsObject * o, FILE * fp)
 {
   GfsInitFraction * init = GFS_INIT_FRACTION (o);
 
-  if (GTS_OBJECT_CLASS (gfs_init_fraction_class ())->parent_class->write)
-    (* GTS_OBJECT_CLASS (gfs_init_fraction_class ())->parent_class->write) 
-      (o, fp);
-  fprintf (fp, " %s { ", init->c->name);
-  gts_surface_write (init->surface, fp);
-  fputs ("}\n", fp);
+
+  (* GTS_OBJECT_CLASS (gfs_init_fraction_class ())->parent_class->write) (o, fp);
+  fprintf (fp, " %s", init->c->name);
+  gfs_surface_write (init->surface, gfs_object_simulation (o), fp);
 }
 
 static gboolean gfs_init_fraction_event (GfsEvent * event, GfsSimulation * sim)
@@ -1728,10 +1681,7 @@ static void gfs_init_fraction_class_init (GfsInitFractionClass * klass)
 
 static void gfs_init_fraction_init (GfsInitFraction * object)
 {
-  object->surface = gts_surface_new (gts_surface_class (),
-				     gts_face_class (),
-				     gts_edge_class (),
-				     gts_vertex_class ());
+  object->surface = GFS_SURFACE (gts_object_new (gfs_surface_class ()));
 }
 
 GfsInitFractionClass * gfs_init_fraction_class (void)
diff --git a/src/event.h b/src/event.h
index 146677a..712b050 100644
--- a/src/event.h
+++ b/src/event.h
@@ -252,7 +252,7 @@ struct _GfsInitFraction {
   GfsGenericInit parent;
 
   GfsVariable * c;
-  GtsSurface * surface;
+  GfsSurface * surface;
 };
 
 typedef struct _GfsInitFractionClass    GfsInitFractionClass;
diff --git a/src/fluid.c b/src/fluid.c
index 66b2df7..9b1b76d 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -2091,140 +2091,6 @@ void gfs_cell_traverse_mixed (FttCell * root,
   cell_traverse_mixed (root, order, flags, func, data);
 }
 
-static void face_overlaps_box (GtsTriangle * t, gpointer * data)
-{
-  GtsBBox * bb = data[0];
-  GtsSurface ** s1 = data[1];
-
-  if (gts_bbox_overlaps_triangle (bb, t)) {
-    if (*s1 == NULL)
-      *s1 = gts_surface_new (gts_surface_class (),
-			     gts_face_class (),
-			     gts_edge_class (),
-			     gts_vertex_class ());
-    gts_surface_add_face (*s1, GTS_FACE (t));
-  }
-}
-
-/**
- * gfs_cell_is_cut:
- * @cell: a #FttCell.
- * @s: a #GtsSurface.
- * @flatten: if set to %TRUE, @cell is flattened in the z direction.
- *
- * Returns: a new #GtsSurface containing the faces of @s which may
- * intersect @cell or %NULL if no faces of @s intersects @cell.
- */
-GtsSurface * gfs_cell_is_cut (FttCell * cell, GtsSurface * s, gboolean flatten)
-{
-  GtsSurface * s1 = NULL;
-  gpointer data[2];
-  GtsBBox bb;
-
-  g_return_val_if_fail (cell != NULL, NULL);
-  g_return_val_if_fail (s != NULL, NULL);
-
-  ftt_cell_bbox (cell, &bb);
-  if (flatten)
-    bb.z1 = bb.z2 = 0.;
-  data[0] = &bb;
-  data[1] = &s1;
-  gts_surface_foreach_face (s, (GtsFunc) face_overlaps_box, data);
-  return s1;
-}
-
-static void cell_traverse_cut (FttCell * cell,
-			       GtsSurface * s,
-			       FttTraverseType order,
-			       FttTraverseFlags flags,
-			       FttCellTraverseCutFunc func,
-			       gpointer data,
-			       gboolean flatten)
-{
-  GtsSurface * s1 = gfs_cell_is_cut (cell, s, flatten && FTT_CELL_IS_LEAF (cell));
-
-  if (s1 == NULL)
-    return;
-  if (order == FTT_PRE_ORDER &&
-      (flags == FTT_TRAVERSE_ALL ||
-       ((flags & FTT_TRAVERSE_LEAFS) != 0 && FTT_CELL_IS_LEAF (cell)) ||
-       ((flags & FTT_TRAVERSE_NON_LEAFS) != 0 && !FTT_CELL_IS_LEAF (cell))))
-    (* func) (cell, s1, data);
-  if (!FTT_CELL_IS_LEAF (cell)) {
-    struct _FttOct * children = cell->children;
-    guint n;
-
-    for (n = 0; n < FTT_CELLS; n++) {
-      FttCell * c = &(children->cell[n]);
-
-      if (!FTT_CELL_IS_DESTROYED (c))
-	cell_traverse_cut (c, s1, order, flags, func, data, flatten);
-    }
-  }
-  if (order == FTT_POST_ORDER &&
-      (flags == FTT_TRAVERSE_ALL ||
-       ((flags & FTT_TRAVERSE_LEAFS) != 0 && FTT_CELL_IS_LEAF (cell)) ||
-       ((flags & FTT_TRAVERSE_NON_LEAFS) != 0 && !FTT_CELL_IS_LEAF (cell))))
-    (* func) (cell, s1, data);
-  gts_object_destroy (GTS_OBJECT (s1));
-}
-
-/**
- * gfs_cell_traverse_cut:
- * @root: the root #FttCell of the tree to traverse.
- * @s: a #GtsSurface.
- * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
- * %FTT_POST_ORDER. 
- * @flags: which types of children are to be visited.
- * @func: the function to call for each visited #FttCell.
- * @data: user data to pass to @func.
- * 
- * Traverses a cell tree starting at the given root #FttCell. Calls
- * the given function for each cell cut by @s.
- */
-void gfs_cell_traverse_cut (FttCell * root,
-			    GtsSurface * s,
-			    FttTraverseType order,
-			    FttTraverseFlags flags,
-			    FttCellTraverseCutFunc func,
-			    gpointer data)
-{
-  g_return_if_fail (root != NULL);
-  g_return_if_fail (s != NULL);
-  g_return_if_fail (func != NULL);
-
-  cell_traverse_cut (root, s, order, flags, func, data, FALSE);
-}
-
-/**
- * gfs_cell_traverse_cut_2D:
- * @root: the root #FttCell of the tree to traverse.
- * @s: a #GtsSurface.
- * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
- * %FTT_POST_ORDER. 
- * @flags: which types of children are to be visited.
- * @func: the function to call for each visited #FttCell.
- * @data: user data to pass to @func.
- * 
- * Traverses a cell tree starting at the given root #FttCell. Calls
- * the given function for each cell cut by @s.
- *
- * The cells are "flattened" in the z-direction.
- */
-void gfs_cell_traverse_cut_2D (FttCell * root,
-			       GtsSurface * s,
-			       FttTraverseType order,
-			       FttTraverseFlags flags,
-			       FttCellTraverseCutFunc func,
-			       gpointer data)
-{
-  g_return_if_fail (root != NULL);
-  g_return_if_fail (s != NULL);
-  g_return_if_fail (func != NULL);
-
-  cell_traverse_cut (root, s, order, flags, func, data, TRUE);
-}
-
 #if FTT_2D
 static FttDirection corner[4][FTT_DIMENSION] = {
   { FTT_LEFT,  FTT_BOTTOM },
diff --git a/src/fluid.h b/src/fluid.h
index 2a77ce3..bc1f09c 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -202,24 +202,6 @@ void                  gfs_cell_traverse_mixed       (FttCell * root,
 						     FttTraverseFlags flags,
 						     FttCellTraverseFunc func,
 						     gpointer data);
-GtsSurface *          gfs_cell_is_cut               (FttCell * cell,
-						     GtsSurface * s,
-						     gboolean flatten);
-typedef void       (* FttCellTraverseCutFunc)       (FttCell * cell,
-						     GtsSurface * s,
-						     gpointer data);
-void                  gfs_cell_traverse_cut         (FttCell * root,
-						     GtsSurface * s,
-						     FttTraverseType order,
-						     FttTraverseFlags flags,
-						     FttCellTraverseCutFunc func,
-						     gpointer data);
-void                  gfs_cell_traverse_cut_2D      (FttCell * root,
-						     GtsSurface * s,
-						     FttTraverseType order,
-						     FttTraverseFlags flags,
-						     FttCellTraverseCutFunc func,
-						     gpointer data);
 gdouble               gfs_interpolate               (FttCell * cell,
 						     FttVector p,
 						     GfsVariable * v);
diff --git a/src/gfs.h b/src/gfs.h
index 4dee384..3bdb2ff 100644
--- a/src/gfs.h
+++ b/src/gfs.h
@@ -38,6 +38,7 @@
 #include <gerris/source.h>
 #include <gerris/vof.h>
 #include <gerris/cartesian.h>
+#include <gerris/surface.h>
 #include <gerris/version.h>
 
 #endif /* GFS_H */
diff --git a/src/init.c b/src/init.c
index 8497f11..09481fc 100644
--- a/src/init.c
+++ b/src/init.c
@@ -134,7 +134,7 @@ GtsObjectClass ** gfs_classes (void)
         gfs_variable_position_class (),
       gfs_variable_distance_class (),
 
-    gfs_surface_class (),
+    gfs_solid_class (),
 
     gfs_init_class (),
     gfs_init_flow_constant_class (),
diff --git a/src/output.c b/src/output.c
index 630b782..b22e73b 100644
--- a/src/output.c
+++ b/src/output.c
@@ -1245,7 +1245,7 @@ static gboolean output_simulation_event (GfsEvent * event, GfsSimulation * sim)
     g_slist_free (domain->variables_io);
     domain->variables_io = output->var;
     domain->binary =       output->binary;
-    sim->output_surface =  output->surface;
+    sim->output_solid   =  output->solid;
     switch (output->format) {
     case GFS:
       gfs_simulation_write (sim,
@@ -1273,7 +1273,7 @@ static gboolean output_simulation_event (GfsEvent * event, GfsSimulation * sim)
     }
     domain->variables_io = NULL;
     domain->binary =       TRUE;
-    sim->output_surface =  TRUE;
+    sim->output_solid   =  TRUE;
     fflush (GFS_OUTPUT (event)->file->fp);
     return TRUE;
   }
@@ -1300,8 +1300,8 @@ static void output_simulation_write (GtsObject * o, FILE * fp)
   }
   if (!output->binary)
     fputs (" binary = 0", fp);
-  if (!output->surface)
-    fputs (" surface = 0", fp);
+  if (!output->solid)
+    fputs (" solid = 0", fp);
   if (output->format == GFS_TEXT)
     fputs (" format = text", fp);
   fputs (" }", fp);
@@ -1309,52 +1309,15 @@ static void output_simulation_write (GtsObject * o, FILE * fp)
 
 static void output_simulation_read (GtsObject ** o, GtsFile * fp)
 {
-  GtsFileVariable var[] = {
-    {GTS_INT,    "depth",    TRUE},
-    {GTS_STRING, "variables",TRUE},
-    {GTS_INT,    "binary",   TRUE},
-    {GTS_INT,    "surface",  TRUE},
-    {GTS_STRING, "format",   TRUE},
-    {GTS_NONE}
-  };
-  gchar * variables = NULL, * format = NULL;
-  GfsOutputSimulation * output = GFS_OUTPUT_SIMULATION (*o);
-  GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (output));
-
-  (* GTS_OBJECT_CLASS (gfs_output_simulation_class ())->parent_class->read) 
-    (o, fp);
+  (* GTS_OBJECT_CLASS (gfs_output_simulation_class ())->parent_class->read) (o, fp);
   if (fp->type == GTS_ERROR)
     return;
 
-  var[0].data = &output->max_depth;
-  var[1].data = &variables;
-  var[2].data = &output->binary;
-  var[3].data = &output->surface;
-  var[4].data = &format;
-  gts_file_assign_variables (fp, var);
-  if (fp->type == GTS_ERROR) {
-    g_free (variables);
-    return;
-  }
-
-  if (variables != NULL) {
-    gchar * error = NULL;
-    GSList * vars = gfs_variables_from_list (domain->variables, variables, &error);
-
-    if (vars == NULL) {
-      gts_file_variable_error (fp, var, "variables",
-			       "unknown variable `%s'", error);
-      g_free (variables);
-      return;
-    }
-    if (output->var)
-      g_slist_free (output->var);
-    output->var = vars;
-    g_free (variables);
-  }
-  else if (output->var == NULL) {
+  GfsOutputSimulation * output = GFS_OUTPUT_SIMULATION (*o);
+  GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (output));
+  if (output->var == NULL) {
     GSList * i = domain->variables;
-
+    
     while (i) {
       if (GFS_VARIABLE1 (i->data)->name)
 	output->var = g_slist_append (output->var, i->data);
@@ -1362,18 +1325,57 @@ static void output_simulation_read (GtsObject ** o, GtsFile * fp)
     }
   }
 
-  if (format != NULL) {
-    if (!strcmp (format, "gfs"))
-      output->format = GFS;
-    else if (!strcmp (format, "text"))
-      output->format = GFS_TEXT;
-    else {
-      gts_file_variable_error (fp, var, "format",
-			       "unknown format `%s'", format);
-      g_free (format);
+  if (fp->type == '{') {
+    GtsFileVariable var[] = {
+      {GTS_INT,    "depth",    TRUE},
+      {GTS_STRING, "variables",TRUE},
+      {GTS_INT,    "binary",   TRUE},
+      {GTS_INT,    "solid",    TRUE},
+      {GTS_STRING, "format",   TRUE},
+      {GTS_NONE}
+    };
+    gchar * variables = NULL, * format = NULL;
+
+    var[0].data = &output->max_depth;
+    var[1].data = &variables;
+    var[2].data = &output->binary;
+    var[3].data = &output->solid;
+    var[4].data = &format;
+    gts_file_assign_variables (fp, var);
+    if (fp->type == GTS_ERROR) {
+      g_free (variables);
       return;
     }
-    g_free (format);
+
+    if (variables != NULL) {
+      gchar * error = NULL;
+      GSList * vars = gfs_variables_from_list (domain->variables, variables, &error);
+
+      if (vars == NULL) {
+	gts_file_variable_error (fp, var, "variables",
+				 "unknown variable `%s'", error);
+	g_free (variables);
+	return;
+      }
+      if (output->var)
+	g_slist_free (output->var);
+      output->var = vars;
+      g_free (variables);
+    }
+
+    if (format != NULL) {
+      if (!strcmp (format, "gfs"))
+	output->format = GFS;
+      else if (!strcmp (format, "text"))
+	output->format = GFS_TEXT;
+      else {
+	gts_file_variable_error (fp, var, "format",
+				 "unknown format `%s'", format);
+	g_free (format);
+	return;
+      }
+      g_free (format);
+    }
   }
 }
 
@@ -1390,7 +1392,7 @@ static void gfs_output_simulation_init (GfsOutputSimulation * object)
   object->max_depth = -1;
   object->var = NULL;
   object->binary = 1;
-  object->surface = 1;
+  object->solid = 1;
   object->format = GFS;
 }
 
diff --git a/src/output.h b/src/output.h
index bc2828b..18b014f 100644
--- a/src/output.h
+++ b/src/output.h
@@ -151,7 +151,7 @@ struct _GfsOutputSimulation {
 
   gint max_depth;
   GSList * var;
-  gboolean binary, surface;
+  gboolean binary, solid;
   GfsOutputSimulationFormat format;
 };
 
diff --git a/src/refine.c b/src/refine.c
index 2b72ee1..f3f61df 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -156,26 +156,34 @@ static void refine_solid_destroy (GtsObject * object)
   (* GTS_OBJECT_CLASS (gfs_refine_solid_class ())->parent_class->destroy) (object);
 }
 
-static void max_kappa (GtsVertex * v, gpointer * data)
+typedef struct {
+  GtsSurface * s;
+  gdouble kappa;
+} KappaData;
+
+static void max_kappa (GtsVertex * v, KappaData * d)
 {
-  GtsSurface * s = data[0];
-  gdouble * max = data[1];
   GtsVector Kh;
 
-  if (gts_vertex_mean_curvature_normal (v, s, Kh)) {
+  if (gts_vertex_mean_curvature_normal (v, d->s, Kh)) {
     gdouble kappa = gts_vector_norm (Kh)/(FTT_DIMENSION - 1);
-    if (kappa > *max)
-      *max = kappa;
+    if (kappa > d->kappa)
+      d->kappa = kappa;
   }
 }
 
 static gdouble solid_curvature (FttCell * cell, FttCellFace * face, 
-				GfsDomain * domain, GtsSurface * s)
-{
-  gdouble kappa = gfs_solid_is_thin (cell, s) ? 1./ftt_cell_size (cell) : 0.;
-  gpointer data[] = {GTS_OBJECT (s)->reserved, &kappa};
-  gts_surface_foreach_vertex (s, (GtsFunc) max_kappa, data);
-  return kappa;
+				GfsDomain * domain, GfsSurface * s)
+{
+  if (s->s) {
+    KappaData d;
+    d.s = s->s;
+    d.kappa = gfs_solid_is_thin (cell, s) ? 1./ftt_cell_size (cell) : 0.;
+    gts_surface_foreach_vertex (d.s, (GtsFunc) max_kappa, &d);
+    return d.kappa;
+  }
+  else /* fixme: need to compute curvature for other types of surfaces */
+    return 0.;
 }
 
 static void refine_solid_read (GtsObject ** o, GtsFile * fp)
@@ -196,10 +204,10 @@ static void refine_solid_read (GtsObject ** o, GtsFile * fp)
 typedef struct {
   GfsRefine * refine;
   GfsDomain * domain;
-  GtsSurface * surface;
+  GfsSurface * surface;
 } RefineCut;
 
-static void refine_cut_cell (FttCell * cell, GtsSurface * s, RefineCut * p)
+static void refine_cut_cell (FttCell * cell, GfsSurface * s, RefineCut * p)
 {
   GTS_OBJECT (s)->reserved = p->surface;
   GFS_REFINE_SOLID (p->refine)->v->data = s;
@@ -210,16 +218,18 @@ static void refine_cut_cell (FttCell * cell, GtsSurface * s, RefineCut * p)
 
 static void gfs_refine_solid_refine (GfsRefine * refine, GfsSimulation * sim)
 {
-  GtsSurface * surface = gfs_simulation_get_surface (sim);
-  if (surface) {
+  if (sim->solids) {
     RefineCut p;
     p.refine = refine;
     p.domain = GFS_DOMAIN (sim);
-    p.surface = surface;
-    gfs_domain_traverse_cut (GFS_DOMAIN (sim), surface, 
-			     FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
-			     (FttCellTraverseCutFunc) refine_cut_cell, &p);
-    gts_object_destroy (GTS_OBJECT (surface));
+    GSList * i = sim->solids->items;
+    while (i) {
+      p.surface = GFS_SOLID (i->data)->s;
+      gfs_domain_traverse_cut (GFS_DOMAIN (sim), p.surface,
+			       FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
+			       (FttCellTraverseCutFunc) refine_cut_cell, &p);
+      i = i->next;
+    }
   }
 }
 
@@ -268,9 +278,7 @@ static void refine_surface_write (GtsObject * o, FILE * fp)
   GfsRefineSurface * d = GFS_REFINE_SURFACE (o);
   
   (* GTS_OBJECT_CLASS (gfs_refine_surface_class ())->parent_class->write) (o, fp);
-  fprintf (fp, " { ");
-  gts_surface_write (d->surface, fp);
-  fputc ('}', fp);
+  gfs_surface_write (d->surface, gfs_object_simulation (o), fp);
   if (d->twod)
     fputs (" { twod = 1 }\n", fp);
   else
@@ -286,45 +294,9 @@ static void refine_surface_read (GtsObject ** o, GtsFile * fp)
     return;
 
   refine = GFS_REFINE_SURFACE (*o);
-  if (fp->type != '{') {
-    FILE * f;
-    GtsFile * gf;
-
-    if (fp->type != GTS_STRING) {
-      gts_file_error (fp, "expecting a string (filename)\n");
-      return;
-    }
-    f = fopen (fp->token->str, "rt");
-    if (f == NULL) {
-      gts_file_error (fp, "cannot open file `%s'\n", fp->token->str);
-      return;
-    }
-    gf = gts_file_new (f);
-    if (gts_surface_read (refine->surface, gf)) {
-      gts_file_error (fp, 
-		      "file `%s' is not a valid GTS file\n"
-		      "%s:%d:%d: %s",
-		      fp->token->str, fp->token->str,
-		      gf->line, gf->pos, gf->error);
-      gts_file_destroy (gf);
-      fclose (f);
-      return;
-    }
-    gts_file_destroy (gf);
-    fclose (f);
-  }
-  else { /* embedded GTS file */
-    fp->scope_max++;
-    gts_file_next_token (fp);
-    if (gts_surface_read (refine->surface, fp))
-      return;
-    if (fp->type != '}') {
-      gts_file_error (fp, "expecting a closing brace");
-      return;
-    }
-    fp->scope_max--;
-  }
-  gts_file_next_token (fp);
+  gfs_surface_read (refine->surface, fp);
+  if (fp->type == GTS_ERROR)
+    return;
 
   if (fp->type == '{') {
     GtsFileVariable var[] = {
@@ -364,10 +336,7 @@ static void gfs_refine_surface_class_init (GfsRefineClass * klass)
 
 static void refine_surface_init (GfsRefineSurface * r)
 {
-  r->surface = gts_surface_new (gts_surface_class (), 
-				gts_face_class (), 
-				gts_edge_class (), 
-				gts_vertex_class ());
+  r->surface = GFS_SURFACE (gts_object_new (gfs_surface_class ()));
 }
 
 GfsRefineClass * gfs_refine_surface_class (void)
@@ -434,8 +403,13 @@ static void refine_distance_read (GtsObject ** o, GtsFile * fp)
   (* GTS_OBJECT_CLASS (gfs_refine_distance_class ())->parent_class->read) (o, fp);
   if (fp->type == GTS_ERROR)
     return;
+  
+  if (!GFS_REFINE_SURFACE (*o)->surface->s) {
+    gts_file_error (fp, "RefineDistance only works with GTS surfaces");
+    return;
+  }
 
-  GFS_REFINE_DISTANCE (*o)->stree = gts_bb_tree_surface (GFS_REFINE_SURFACE (*o)->surface);
+  GFS_REFINE_DISTANCE (*o)->stree = gts_bb_tree_surface (GFS_REFINE_SURFACE (*o)->surface->s);
 }
 
 static void gfs_refine_distance_class_init (GfsRefineClass * klass)
@@ -498,7 +472,7 @@ static gdouble cell_height (FttCell * cell,
 {
   FttVector pos;
   ftt_cell_pos (cell, &pos);
-  return interpolated_value (refine->surface, &pos);
+  return interpolated_value (refine->surface->s, &pos);
 }
 
 static void refine_height_read (GtsObject ** o, GtsFile * fp)
@@ -513,6 +487,13 @@ static void refine_height_read (GtsObject ** o, GtsFile * fp)
   }
 
   (* GTS_OBJECT_CLASS (gfs_refine_height_class ())->parent_class->read) (o, fp);
+  if (fp->type == GTS_ERROR)
+    return;
+
+  if (!GFS_REFINE_SURFACE (*o)->surface->s) {
+    gts_file_error (fp, "RefineDistance only works with GTS surfaces");
+    return;
+  }
 }
 
 static void gfs_refine_height_class_init (GfsRefineClass * klass)
diff --git a/src/refine.h b/src/refine.h
index e63646f..4a4e6e8 100644
--- a/src/refine.h
+++ b/src/refine.h
@@ -75,7 +75,7 @@ typedef struct _GfsRefineSurface         GfsRefineSurface;
 struct _GfsRefineSurface {
   GfsRefine parent;
 
-  GtsSurface * surface;
+  GfsSurface * surface;
   gboolean twod;
 };
 
diff --git a/src/simulation.c b/src/simulation.c
index 4829025..cecc525 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -46,7 +46,7 @@ static void simulation_destroy (GtsObject * object)
 			 (GtsFunc) gts_object_destroy, NULL);
   gts_object_destroy (GTS_OBJECT (sim->events));
   gts_object_destroy (GTS_OBJECT (sim->adapts));
-  gts_object_destroy (GTS_OBJECT (sim->surfaces));
+  gts_object_destroy (GTS_OBJECT (sim->solids));
 
   g_slist_foreach (sim->modules, (GFunc) g_module_close, NULL);
   g_slist_free (sim->modules);
@@ -275,8 +275,8 @@ static void simulation_read (GtsObject ** object, GtsFile * fp)
 	gts_container_add (GTS_CONTAINER (sim->adapts), GTS_CONTAINEE (object));
 	gts_container_add (GTS_CONTAINER (sim->events), GTS_CONTAINEE (object));
       }
-      else if (GFS_IS_SURFACE (object)) {
-	gts_container_add (GTS_CONTAINER (sim->surfaces), GTS_CONTAINEE (object));
+      else if (GFS_IS_SOLID (object)) {
+	gts_container_add (GTS_CONTAINER (sim->solids), GTS_CONTAINEE (object));
 	gts_container_add (GTS_CONTAINER (sim->events), GTS_CONTAINEE (object));
       }
       else if (GFS_IS_EVENT (object))
@@ -298,7 +298,7 @@ static void simulation_read (GtsObject ** object, GtsFile * fp)
   sim->refines->items = g_slist_reverse (sim->refines->items);
   sim->adapts->items = g_slist_reverse (sim->adapts->items);
   sim->events->items = g_slist_reverse (sim->events->items);
-  sim->surfaces->items = g_slist_reverse (sim->surfaces->items);
+  sim->solids->items = g_slist_reverse (sim->solids->items);
   sim->modules = g_slist_reverse (sim->modules);
 }
 
@@ -684,10 +684,10 @@ static void simulation_init (GfsSimulation * object)
   gfs_multilevel_params_init (&object->projection_params);
   gfs_multilevel_params_init (&object->approx_projection_params);
 
-  object->surfaces = GTS_SLIST_CONTAINER (gts_container_new
-					  (GTS_CONTAINER_CLASS
-					   (gts_slist_container_class ())));
-  object->output_surface = TRUE;
+  object->solids = GTS_SLIST_CONTAINER (gts_container_new
+					(GTS_CONTAINER_CLASS
+					 (gts_slist_container_class ())));
+  object->output_solid = TRUE;
 
   object->refines = GTS_SLIST_CONTAINER (gts_container_new
 					 (GTS_CONTAINER_CLASS
@@ -743,34 +743,6 @@ static void init_non_variable (GfsEvent * event, GfsSimulation * sim)
 }
 
 /**
- * gfs_simulation_get_surface:
- * @sim: a #GfsSimulation.
- *
- * Returns: a new #GtsSurface containing all the solid surfaces
- * defined in @sim, or %NULL if @sim does not contain any solid
- * surface.
- */
-GtsSurface * gfs_simulation_get_surface (GfsSimulation * sim)
-{
-  g_return_val_if_fail (sim != NULL, NULL);
-
-  if (sim->surfaces->items == NULL)
-    return NULL;
-  else {
-    GtsSurface * s = gts_surface_new (gts_surface_class (), 
-				      gts_face_class (), 
-				      gts_edge_class (), 
-				      gts_vertex_class ());
-    GSList * i = sim->surfaces->items;
-    while (i) {
-      gts_surface_merge (s, GFS_SURFACE (i->data)->s);
-      i = i->next;
-    }
-    return s;
-  }
-}
-
-/**
  * gfs_simulation_init:
  * @sim: a #GfsSimulation.
  *
@@ -832,6 +804,25 @@ static void set_permanent (FttCell * cell)
 }
 
 /**
+ * gfs_simulation_get_solids:
+ * @sim: a #GfsSimulation.
+ *
+ * Returns: a new list of #GfsSurface defining the solid boundaries
+ * contained in @sim.
+ */
+GSList * gfs_simulation_get_solids (GfsSimulation * sim)
+{
+  g_return_val_if_fail (sim != NULL, NULL);
+
+  GSList * solids = NULL, * i = sim->solids->items;
+  while (i) {
+    solids = g_slist_prepend (solids, GFS_SOLID (i->data)->s);
+    i = i->next;
+  }
+  return solids;
+}
+
+/**
  * gfs_simulation_refine:
  * @sim: a #GfsSimulation.
  *
@@ -869,13 +860,13 @@ void gfs_simulation_refine (GfsSimulation * sim)
   gfs_domain_match (domain);
   gfs_domain_timer_stop (domain, "simulation_refine");
 
-  GtsSurface * surface = gfs_simulation_get_surface (sim);
-  if (surface) {
+  GSList * solids = gfs_simulation_get_solids (sim);
+  if (solids) {
     gfs_domain_timer_start (domain, "solid_fractions");
-    sim->thin = gfs_domain_init_solid_fractions (domain, surface, TRUE,
+    sim->thin = gfs_domain_init_solid_fractions (domain, solids, TRUE,
 						 (FttCellCleanupFunc) gfs_cell_cleanup, NULL,  
 						 NULL);
-    gts_object_destroy (GTS_OBJECT (surface));
+    g_slist_free (solids);
     gfs_domain_match (domain);
     gfs_domain_traverse_mixed (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
 			       (FttCellTraverseFunc) set_permanent, NULL);
diff --git a/src/simulation.h b/src/simulation.h
index bf9a3a1..870608e 100644
--- a/src/simulation.h
+++ b/src/simulation.h
@@ -71,9 +71,9 @@ struct _GfsSimulation {
   GtsSListContainer * events;
   GSList * modules, * variables, * globals;
 
-  GtsSListContainer * surfaces;
+  GtsSListContainer * solids;
   guint thin;
-  gboolean output_surface;
+  gboolean output_solid;
 
   gdouble tnext;
 };
@@ -95,12 +95,12 @@ struct _GfsSimulationClass {
 
 GfsSimulationClass * gfs_simulation_class        (void);
 GfsSimulation *      gfs_simulation_new          (GfsSimulationClass * klass);
-GtsSurface *         gfs_simulation_get_surface  (GfsSimulation * sim);
 void                 gfs_simulation_init         (GfsSimulation * sim);
 void                 gfs_simulation_write        (GfsSimulation * sim,
 						  gint max_depth,  
 						  FILE * fp);
 GfsSimulation *      gfs_simulation_read         (GtsFile * fp);
+GSList *             gfs_simulation_get_solids   (GfsSimulation * sim);
 void                 gfs_simulation_refine       (GfsSimulation * sim);
 void                 gfs_simulation_set_timestep (GfsSimulation * sim);
 void                 gfs_time_init               (GfsTime * t);
diff --git a/src/solid.c b/src/solid.c
index 2e637fa..85e2f25 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -221,7 +221,7 @@ static void triangle_face_intersection (GtsTriangle * t, CellFace * f)
   }
 }
 
-static void face_new (CellFace * f, FttCell * cell, GtsSurface * s, FttVector * h)
+static void face_new (CellFace * f, FttCell * cell, GfsSurface * s, FttVector * h)
 {
   FttVector p;
 
@@ -234,7 +234,7 @@ static void face_new (CellFace * f, FttCell * cell, GtsSurface * s, FttVector *
   f->n[0] = f->n[1] = f->n[2] = f->n[3] = 0;
   f->inside[0] = f->inside[1] = f->inside[2] = f->inside[3] = 0;
 
-  gts_surface_foreach_face (s, (GtsFunc) triangle_face_intersection, f);
+  gts_surface_foreach_face (s->s, (GtsFunc) triangle_face_intersection, f);
 }
 
 static gboolean solid_face_is_thin (CellFace * f)
@@ -260,14 +260,14 @@ static gboolean solid_face_is_thin (CellFace * f)
 /**
  * gfs_set_2D_solid_fractions_from_surface:
  * @cell: a #FttCell.
- * @s: a #GtsSurface.
+ * @s: a #GfsSurface.
  *
  * Sets the 2D volume fractions of @cell cut by @s.
  *
  * Returns: %TRUE if the cell is thin, %FALSE otherwise;
  */
 gboolean gfs_set_2D_solid_fractions_from_surface (FttCell * cell,
-						  GtsSurface * s)
+						  GfsSurface * s)
 {
   GfsSolidVector * solid;
   FttVector h;
@@ -353,7 +353,7 @@ static void deal_with_thin_cell (FttCell * cell, InitSolidParams * p)
 #if FTT_2D /* 2D */
 
 static void set_solid_fractions_from_surface (FttCell * cell,
-					      GtsSurface * s,
+					      GfsSurface * s,
 					      InitSolidParams * p)
 {
   if (gfs_set_2D_solid_fractions_from_surface (cell, s)) {
@@ -367,7 +367,7 @@ static void set_solid_fractions_from_surface (FttCell * cell,
 /**
  * gfs_solid_is_thin:
  * @cell: a #FttCell.
- * @s: a #GtsSurface.
+ * @s: a #GfsSurface.
  *
  * @s is "thin" relative to @cell if the miminum distance between
  * non-connected faces of @s cutting @cell is smaller than the size of
@@ -375,7 +375,7 @@ static void set_solid_fractions_from_surface (FttCell * cell,
  *
  * Returns: %TRUE if @s is a thin surface, %FALSE otherwise.
  */
-gboolean gfs_solid_is_thin (FttCell * cell, GtsSurface * s)
+gboolean gfs_solid_is_thin (FttCell * cell, GfsSurface * s)
 {
   CellFace f;
   FttVector h;
@@ -482,7 +482,7 @@ static guint topology (CellCube * cube)
   return nl;
 }
 
-static void cube_new (CellCube * cube, FttCell * cell, GtsSurface * s, FttVector * o, FttVector * h)
+static void cube_new (CellCube * cube, FttCell * cell, GfsSurface * s, FttVector * o, FttVector * h)
 {
   guint i;
 
@@ -497,10 +497,12 @@ static void cube_new (CellCube * cube, FttCell * cell, GtsSurface * s, FttVector
     cube->p[i].z = o->z + h->z*vertex[i].z;
   }
 
-  gts_surface_foreach_face (s, (GtsFunc) triangle_cube_intersection, cube);  
+  gts_surface_foreach_face (s->s, (GtsFunc) triangle_cube_intersection, cube);  
 }
 
-static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s, InitSolidParams * p)
+static void set_solid_fractions_from_surface (FttCell * cell, 
+					      GfsSurface * s, 
+					      InitSolidParams * p)
 {
   GfsSolidVector * solid = GFS_STATE (cell)->solid;
   CellCube cube;
@@ -647,7 +649,7 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s, In
 /**
  * gfs_solid_is_thin:
  * @cell: a #FttCell.
- * @s: a #GtsSurface.
+ * @s: a #GfsSurface.
  *
  * @s is "thin" relative to @cell if the miminum distance between
  * non-connected faces of @s cutting @cell is smaller than the size of
@@ -655,7 +657,7 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s, In
  *
  * Returns: %TRUE if @s is a thin surface, %FALSE otherwise.
  */
-gboolean gfs_solid_is_thin (FttCell * cell, GtsSurface * s)
+gboolean gfs_solid_is_thin (FttCell * cell, GfsSurface * s)
 {
   CellCube cube;
   FttVector o, h;
@@ -931,7 +933,7 @@ static void match_fractions (FttCell * cell, GfsVariable * status)
 /**
  * gfs_domain_init_solid_fractions:
  * @domain: a #GfsDomain.
- * @s: an orientable surface defining the solid boundary.
+ * @i: a list of #GfsSurfaces.
  * @destroy_solid: controls what to do with solid cells.
  * @cleanup: a #FttCellCleanupFunc or %NULL.
  * @data: user data to pass to @cleanup.
@@ -945,7 +947,7 @@ static void match_fractions (FttCell * cell, GfsVariable * status)
  * Returns: the number of thin cells.
  */
 guint gfs_domain_init_solid_fractions (GfsDomain * domain,
-				       GtsSurface * s,
+				       GSList * i,
 				       gboolean destroy_solid,
 				       FttCellCleanupFunc cleanup,
 				       gpointer data,
@@ -954,7 +956,6 @@ guint gfs_domain_init_solid_fractions (GfsDomain * domain,
   InitSolidParams p;
 
   g_return_val_if_fail (domain != NULL, 0);
-  g_return_val_if_fail (s != NULL, 0);
 
   p.destroy_solid = destroy_solid;
   p.cleanup = cleanup;
@@ -963,8 +964,11 @@ guint gfs_domain_init_solid_fractions (GfsDomain * domain,
   p.thin = 0;
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 			    (FttCellTraverseFunc) gfs_cell_reset, p.status);
-  gfs_domain_traverse_cut (domain, s, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
-			   (FttCellTraverseCutFunc) set_solid_fractions_from_surface, &p);
+  while (i) {
+    gfs_domain_traverse_cut (domain, i->data, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS,
+			     (FttCellTraverseCutFunc) set_solid_fractions_from_surface, &p);
+    i = i->next;
+  }
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttCellTraverseFunc) paint_mixed_leaf, p.status);
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
@@ -1168,14 +1172,14 @@ static void restore_solid (FttCell * cell, gpointer * data)
 /**
  * gfs_domain_init_fraction:
  * @domain: a #GfsDomain.
- * @s: an orientable surface defining the interface boundary.
+ * @s: a surface defining the interface boundary.
  * @c: a #GfsVariable.
  *
  * Initializes the fraction @c of the interface @s contained in all
  * the cells of @domain.
  */
 void gfs_domain_init_fraction (GfsDomain * domain,
-			       GtsSurface * s,
+			       GfsSurface * s,
 			       GfsVariable * c)
 {
   gboolean not_cut = TRUE;
@@ -1190,7 +1194,9 @@ void gfs_domain_init_fraction (GfsDomain * domain,
 
   gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
 			    (FttCellTraverseFunc) save_solid, c);
-  gfs_domain_init_solid_fractions (domain, s, FALSE, NULL, NULL, status);
+  GSList * l = g_slist_prepend (NULL, s);
+  gfs_domain_init_solid_fractions (domain, l, FALSE, NULL, NULL, status);
+  g_slist_free (l);
   data[0] = c;
   data[1] = &not_cut;
   data[2] = status;
@@ -1461,205 +1467,62 @@ void gfs_solid_coarse_fine (FttCell * parent)
 #endif /* 3D */
 }
 
-/* GfsSurface: Object */
-
-static void check_solid_surface (GtsSurface * s, 
-				 const gchar * fname,
-				 GtsFile * fp)
-{
-  GString * name = g_string_new ("surface");
-
-  if (fname) {
-    g_string_append (name, " `");
-    g_string_append (name, fname);
-    g_string_append_c (name, '\'');
-  }
-
-  if (!gts_surface_is_orientable (s))
-    gts_file_error (fp, "%s is not orientable", name->str);
-  g_string_free (name, TRUE);
-}
+/* GfsSolid: Object */
 
-static void gfs_surface_read (GtsObject ** o, GtsFile * fp)
+static void gfs_solid_read (GtsObject ** o, GtsFile * fp)
 {
-  (* GTS_OBJECT_CLASS (gfs_surface_class ())->parent_class->read) (o, fp);
+  (* GTS_OBJECT_CLASS (gfs_solid_class ())->parent_class->read) (o, fp);
   if (fp->type == GTS_ERROR)
     return;
 
-  GtsSurface * surface = GFS_SURFACE (*o)->s;
-  if (fp->type == '{') {
-    fp->scope_max++;
-    gts_file_next_token (fp);
-    if (gts_surface_read (surface, fp))
-      return;
-    if (fp->type != '}') {
-      gts_file_error (fp, "expecting a closing brace");
-      return;
-    }
-    check_solid_surface (surface, NULL, fp);
-    if (fp->type == GTS_ERROR)
-      return;
-    fp->scope_max--;
-  }
-  else {
-    if (fp->type != GTS_STRING) {
-      gts_file_error (fp, "expecting a string (filename)");
-      return;
-    }
-    FILE * fptr = fopen (fp->token->str, "rt");
-    if (fptr == NULL) {
-      gts_file_error (fp, "cannot open file `%s'", fp->token->str);
-      return;
-    }
-    GtsFile * fp1 = gts_file_new (fptr);
-    if (gts_surface_read (surface, fp1)) {
-      gts_file_error (fp, 
-		      "file `%s' is not a valid GTS file\n"
-		      "%s:%d:%d: %s",
-		      fp->token->str, fp->token->str,
-		      fp1->line, fp1->pos, fp1->error);
-      gts_file_destroy (fp1);
-      fclose (fptr);
-      return;
-    }
-    gts_file_destroy (fp1);
-    fclose (fptr);
-  
-    check_solid_surface (surface, fp->token->str, fp);
-    if (fp->type == GTS_ERROR)
-      return;
-  }
-  gts_file_next_token (fp);
-
-  if (fp->type == '{') {
-    GtsVector r = {0.,0.,0.}, s = {1.,1.,1.}, t = {0.,0.,0.};
-    gdouble angle = 0., scale = 1.;
-    gboolean flip = FALSE;
-    GtsFileVariable var[] = {
-      {GTS_DOUBLE, "rx", TRUE},
-      {GTS_DOUBLE, "ry", TRUE},
-      {GTS_DOUBLE, "rz", TRUE},
-      {GTS_DOUBLE, "sx", TRUE},
-      {GTS_DOUBLE, "sy", TRUE},
-      {GTS_DOUBLE, "sz", TRUE},
-      {GTS_DOUBLE, "tx", TRUE},
-      {GTS_DOUBLE, "ty", TRUE},
-      {GTS_DOUBLE, "tz", TRUE},
-      {GTS_DOUBLE, "scale", TRUE},
-      {GTS_DOUBLE, "angle", TRUE},
-      {GTS_INT,  "flip", TRUE},
-      {GTS_NONE}
-    };
-    GtsFileVariable * v = var;
-
-    (v++)->data = &r[0];
-    (v++)->data = &r[1];
-    (v++)->data = &r[2];
-
-    (v++)->data = &s[0];
-    (v++)->data = &s[1];
-    (v++)->data = &s[2];
-
-    (v++)->data = &t[0];
-    (v++)->data = &t[1];
-    (v++)->data = &t[2];
-
-    (v++)->data = &scale;
-    (v++)->data = &angle;
-
-    (v++)->data = &flip;
-
-    gts_file_assign_variables (fp, var);
-    if (fp->type == GTS_ERROR)
-      return;
-
-    if (var[9].set)
-      s[0] = s[1] = s[2] = scale;
-    if (var[10].set && gts_vector_norm (r) == 0.) {
-      gts_file_variable_error (fp, var, "angle",
-			       "a non-zero rotation vector must be specified");
-      return;
-    }
-    
-    GtsMatrix * m = gts_matrix_translate (NULL, t);
-    if (angle != 0.) {
-      GtsMatrix * mr = gts_matrix_rotate (NULL, r, angle*M_PI/180.);
-      GtsMatrix * m1 = gts_matrix_product (m, mr);
-      gts_matrix_destroy (m);
-      gts_matrix_destroy (mr);
-      m = m1;
-    }
-    GtsMatrix * ms = gts_matrix_scale (NULL, s);
-    GtsMatrix * M = gts_matrix_product (m, ms);
-    gts_matrix_destroy (m);
-    gts_matrix_destroy (ms);
-    gts_surface_foreach_vertex (surface, (GtsFunc) gts_point_transform, M);
-    gts_matrix_destroy (M);
-
-    if (flip)
-      gts_surface_foreach_face (surface, (GtsFunc) gts_triangle_revert, NULL);
-  }
+  gfs_surface_read (GFS_SOLID (*o)->s, fp);
 }
 
-static void gfs_surface_write (GtsObject * o, FILE * fp)
+static void gfs_solid_write (GtsObject * o, FILE * fp)
 {
   GfsSimulation * sim = gfs_object_simulation (o);
-  if (sim->output_surface) {
-    (* GTS_OBJECT_CLASS (gfs_surface_class ())->parent_class->write) (o, fp);
-
-    fputs (" { ", fp);
-    GtsSurface * s = GFS_SURFACE (o)->s;
-    if (GFS_DOMAIN (sim)->binary) {
-      gboolean binary = GTS_POINT_CLASS (s->vertex_class)->binary;
-      GTS_POINT_CLASS (s->vertex_class)->binary = TRUE;
-      gts_surface_write (s, fp);
-      GTS_POINT_CLASS (s->vertex_class)->binary = binary;
-    }
-    else
-      gts_surface_write (s, fp);
-    fputc ('}', fp);
+  if (sim->output_solid) {
+    (* GTS_OBJECT_CLASS (gfs_solid_class ())->parent_class->write) (o, fp);
+    gfs_surface_write (GFS_SOLID (o)->s, sim, fp);
   }
 }
 
-static void gfs_surface_destroy (GtsObject * object)
+static void gfs_solid_destroy (GtsObject * object)
 {
-  gts_object_destroy (GTS_OBJECT (GFS_SURFACE (object)->s));
+  gts_object_destroy (GTS_OBJECT (GFS_SOLID (object)->s));
 
-  (* GTS_OBJECT_CLASS (gfs_surface_class ())->parent_class->destroy) (object);
+  (* GTS_OBJECT_CLASS (gfs_solid_class ())->parent_class->destroy) (object);
 }
 
-static void gfs_surface_class_init (GtsSurfaceClass * klass)
+static void gfs_solid_class_init (GtsObjectClass * klass)
 {
-  GTS_OBJECT_CLASS (klass)->read = gfs_surface_read;
-  GTS_OBJECT_CLASS (klass)->write = gfs_surface_write;
-  GTS_OBJECT_CLASS (klass)->destroy = gfs_surface_destroy;
+  klass->read = gfs_solid_read;
+  klass->write = gfs_solid_write;
+  klass->destroy = gfs_solid_destroy;
 }
 
-static void gfs_surface_init (GfsSurface * object)
+static void gfs_solid_init (GfsSolid * object)
 {
-  object->s = gts_surface_new (gts_surface_class (), 
-			       gts_face_class (), 
-			       gts_edge_class (), 
-			       gts_vertex_class ());
+  object->s = GFS_SURFACE (gts_object_new (gfs_surface_class ()));
   GFS_EVENT (object)->istep = G_MAXINT/2;
 }
 
-GfsEventClass * gfs_surface_class (void)
+GfsEventClass * gfs_solid_class (void)
 {
   static GfsEventClass * klass = NULL;
 
   if (klass == NULL) {
-    GtsObjectClassInfo gfs_surface_info = {
-      "GfsSurface",
-      sizeof (GfsSurface),
+    GtsObjectClassInfo gfs_solid_info = {
+      "GfsSolid",
+      sizeof (GfsSolid),
       sizeof (GfsEventClass),
-      (GtsObjectClassInitFunc) gfs_surface_class_init,
-      (GtsObjectInitFunc) gfs_surface_init,
+      (GtsObjectClassInitFunc) gfs_solid_class_init,
+      (GtsObjectInitFunc) gfs_solid_init,
       (GtsArgSetFunc) NULL,
       (GtsArgGetFunc) NULL
     };
     klass = gts_object_class_new (GTS_OBJECT_CLASS (gfs_event_class ()),
-				  &gfs_surface_info);
+				  &gfs_solid_info);
   }
 
   return klass;
diff --git a/src/solid.h b/src/solid.h
index a9c6d61..bb992e6 100644
--- a/src/solid.h
+++ b/src/solid.h
@@ -31,11 +31,11 @@ extern "C" {
 
 void         gfs_cell_fluid                              (FttCell * cell);
 gboolean     gfs_solid_is_thin                           (FttCell * cell, 
-							  GtsSurface * s);
+							  GfsSurface * s);
 gboolean     gfs_set_2D_solid_fractions_from_surface     (FttCell * cell,
-							  GtsSurface * s);
+							  GfsSurface * s);
 guint        gfs_domain_init_solid_fractions             (GfsDomain * domain,
-							  GtsSurface * s,
+							  GSList * i,
 							  gboolean destroy_solid,
 							  FttCellCleanupFunc cleanup,
 							  gpointer data,
@@ -43,7 +43,7 @@ guint        gfs_domain_init_solid_fractions             (GfsDomain * domain,
 void         gfs_cell_init_solid_fractions_from_children (FttCell * cell);
 gboolean     gfs_cell_check_solid_fractions              (FttCell * root);
 void         gfs_domain_init_fraction                    (GfsDomain * domain,
-							  GtsSurface * s,
+							  GfsSurface * s,
 							  GfsVariable * c);
 void         gfs_cell_cm                                 (const FttCell * cell, 
 							  FttVector * cm);
@@ -53,25 +53,25 @@ void         gfs_face_ca                                 (const FttCellFace * fa
 							  FttVector * ca);
 void         gfs_solid_coarse_fine                       (FttCell * parent);
 
-/* GfsSurface: Header */
+/* GfsSolid: Header */
 
-typedef struct _GfsSurface         GfsSurface;
+typedef struct _GfsSolid         GfsSolid;
 
-struct _GfsSurface {
+struct _GfsSolid {
   /*< private >*/
   GfsEvent parent;
 
   /*< public >*/
-  GtsSurface * s;
+  GfsSurface * s;
 };
 
-#define GFS_SURFACE(obj)            GTS_OBJECT_CAST (obj,\
-					         GfsSurface,\
-					         gfs_surface_class ())
-#define GFS_IS_SURFACE(obj)         (gts_object_is_from_class (obj,\
-						 gfs_surface_class ()))
+#define GFS_SOLID(obj)            GTS_OBJECT_CAST (obj,\
+					         GfsSolid,\
+					         gfs_solid_class ())
+#define GFS_IS_SOLID(obj)         (gts_object_is_from_class (obj,\
+						 gfs_solid_class ()))
 
-GfsEventClass * gfs_surface_class  (void);
+GfsEventClass * gfs_solid_class  (void);
 
 #ifdef __cplusplus
 }
diff --git a/src/surface.c b/src/surface.c
new file mode 100644
index 0000000..9a5605f
--- /dev/null
+++ b/src/surface.c
@@ -0,0 +1,392 @@
+/* Gerris - The GNU Flow Solver
+ * Copyright (C) 2001 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 "surface.h"
+#include "simulation.h"
+
+/* GfsSurface: Object */
+
+static void check_solid_surface (GtsSurface * s, 
+				 const gchar * fname,
+				 GtsFile * fp)
+{
+  GString * name = g_string_new ("surface");
+
+  if (fname) {
+    g_string_append (name, " `");
+    g_string_append (name, fname);
+    g_string_append_c (name, '\'');
+  }
+
+  if (!gts_surface_is_orientable (s))
+    gts_file_error (fp, "%s is not orientable", name->str);
+  g_string_free (name, TRUE);
+}
+
+static void surface_read (GtsObject ** o, GtsFile * fp)
+{
+  GfsSurface * surface = GFS_SURFACE (*o);
+
+  if (fp->type == '{') {
+    fp->scope_max++;
+    gts_file_next_token (fp);
+    g_assert (!surface->s);
+    surface->s = gts_surface_new (gts_surface_class (), 
+				  gts_face_class (), 
+				  gts_edge_class (), 
+				  gts_vertex_class ());
+    if (gts_surface_read (surface->s, fp))
+      return;
+    if (fp->type != '}') {
+      gts_file_error (fp, "expecting a closing brace");
+      return;
+    }
+    check_solid_surface (surface->s, NULL, fp);
+    if (fp->type == GTS_ERROR)
+      return;
+    fp->scope_max--;
+  }
+  else {
+    if (fp->type != GTS_STRING) {
+      gts_file_error (fp, "expecting a string (filename)");
+      return;
+    }
+    FILE * fptr = fopen (fp->token->str, "rt");
+    if (fptr == NULL) {
+      gts_file_error (fp, "cannot open file `%s'", fp->token->str);
+      return;
+    }
+    GtsFile * fp1 = gts_file_new (fptr);
+    surface->s = gts_surface_new (gts_surface_class (), 
+				  gts_face_class (), 
+				  gts_edge_class (), 
+				  gts_vertex_class ());
+    if (gts_surface_read (surface->s, fp1)) {
+      gts_file_error (fp, 
+		      "file `%s' is not a valid GTS file\n"
+		      "%s:%d:%d: %s",
+		      fp->token->str, fp->token->str,
+		      fp1->line, fp1->pos, fp1->error);
+      gts_file_destroy (fp1);
+      fclose (fptr);
+      return;
+    }
+    gts_file_destroy (fp1);
+    fclose (fptr);
+    check_solid_surface (surface->s, fp->token->str, fp);
+    if (fp->type == GTS_ERROR)
+      return;
+  }
+  gts_file_next_token (fp);
+
+  if (fp->type == '{') {
+    GtsVector r = {0.,0.,0.}, s = {1.,1.,1.}, t = {0.,0.,0.};
+    gdouble angle = 0., scale = 1.;
+    gboolean flip = FALSE;
+    GtsFileVariable var[] = {
+      {GTS_DOUBLE, "rx", TRUE},
+      {GTS_DOUBLE, "ry", TRUE},
+      {GTS_DOUBLE, "rz", TRUE},
+      {GTS_DOUBLE, "sx", TRUE},
+      {GTS_DOUBLE, "sy", TRUE},
+      {GTS_DOUBLE, "sz", TRUE},
+      {GTS_DOUBLE, "tx", TRUE},
+      {GTS_DOUBLE, "ty", TRUE},
+      {GTS_DOUBLE, "tz", TRUE},
+      {GTS_DOUBLE, "scale", TRUE},
+      {GTS_DOUBLE, "angle", TRUE},
+      {GTS_INT,  "flip", TRUE},
+      {GTS_NONE}
+    };
+    GtsFileVariable * v = var;
+
+    (v++)->data = &r[0];
+    (v++)->data = &r[1];
+    (v++)->data = &r[2];
+
+    (v++)->data = &s[0];
+    (v++)->data = &s[1];
+    (v++)->data = &s[2];
+
+    (v++)->data = &t[0];
+    (v++)->data = &t[1];
+    (v++)->data = &t[2];
+
+    (v++)->data = &scale;
+    (v++)->data = &angle;
+
+    (v++)->data = &flip;
+
+    gts_file_assign_variables (fp, var);
+    if (fp->type == GTS_ERROR)
+      return;
+
+    if (var[9].set)
+      s[0] = s[1] = s[2] = scale;
+    if (var[10].set && gts_vector_norm (r) == 0.) {
+      gts_file_variable_error (fp, var, "angle",
+			       "a non-zero rotation vector must be specified");
+      return;
+    }
+    
+    GtsMatrix * m = gts_matrix_translate (NULL, t);
+    if (angle != 0.) {
+      GtsMatrix * mr = gts_matrix_rotate (NULL, r, angle*M_PI/180.);
+      GtsMatrix * m1 = gts_matrix_product (m, mr);
+      gts_matrix_destroy (m);
+      gts_matrix_destroy (mr);
+      m = m1;
+    }
+    GtsMatrix * ms = gts_matrix_scale (NULL, s);
+    GtsMatrix * M = gts_matrix_product (m, ms);
+    gts_matrix_destroy (m);
+    gts_matrix_destroy (ms);
+    gts_surface_foreach_vertex (surface->s, (GtsFunc) gts_point_transform, M);
+    gts_matrix_destroy (M);
+
+    if (flip)
+      gts_surface_foreach_face (surface->s, (GtsFunc) gts_triangle_revert, NULL);
+  }
+}
+
+static void surface_write (GtsObject * o, FILE * fp)
+{
+  fputs (" { ", fp);
+  GtsSurface * s = GFS_SURFACE (o)->s;
+  if (GFS_DOMAIN (gfs_object_simulation (o))->binary) {
+    gboolean binary = GTS_POINT_CLASS (s->vertex_class)->binary;
+    GTS_POINT_CLASS (s->vertex_class)->binary = TRUE;
+    gts_surface_write (s, fp);
+    GTS_POINT_CLASS (s->vertex_class)->binary = binary;
+  }
+  else
+    gts_surface_write (s, fp);
+  fputc ('}', fp);
+}
+
+static void surface_destroy (GtsObject * object)
+{
+  if (GFS_SURFACE (object)->s)
+    gts_object_destroy (GTS_OBJECT (GFS_SURFACE (object)->s));
+
+  (* GTS_OBJECT_CLASS (gfs_surface_class ())->parent_class->destroy) (object);
+}
+
+static void gfs_surface_class_init (GtsObjectClass * klass)
+{
+  klass->read = surface_read;
+  klass->write = surface_write;
+  klass->destroy = surface_destroy;
+}
+
+GtsObjectClass * gfs_surface_class (void)
+{
+  static GtsObjectClass * klass = NULL;
+
+  if (klass == NULL) {
+    GtsObjectClassInfo gfs_surface_info = {
+      "GfsSurface",
+      sizeof (GfsSurface),
+      sizeof (GtsObjectClass),
+      (GtsObjectClassInitFunc) gfs_surface_class_init,
+      (GtsObjectInitFunc) NULL,
+      (GtsArgSetFunc) NULL,
+      (GtsArgGetFunc) NULL
+    };
+    klass = gts_object_class_new (gts_object_class (), &gfs_surface_info);
+  }
+
+  return klass;
+}
+
+/**
+ * gfs_surface_read:
+ * @s: a #GfsSurface.
+ * @fp: a #GtsFile.
+ * 
+ * Calls the read() method of @s.
+ */
+void gfs_surface_read (GfsSurface * s, GtsFile * fp)
+{
+  GtsObject * o = (GtsObject *) s;
+
+  g_return_if_fail (s != NULL);
+  g_return_if_fail (fp != NULL);
+
+  (* GTS_OBJECT (s)->klass->read) (&o, fp);
+}
+
+/**
+ * gfs_surface_write:
+ * @s: a #GfsSurface.
+ * @sim: a #GfsSimulation.
+ * @fp: a file pointer.
+ * 
+ * Calls the write() method of @s.
+ */
+void gfs_surface_write (GfsSurface * s, gpointer sim, FILE * fp)
+{
+  g_return_if_fail (s != NULL);
+  g_return_if_fail (fp != NULL);
+
+  GTS_OBJECT (s)->reserved = sim;
+  (* GTS_OBJECT (s)->klass->write) (GTS_OBJECT (s), fp);
+}
+
+static void face_overlaps_box (GtsTriangle * t, gpointer * data)
+{
+  GtsBBox * bb = data[0];
+  GtsSurface ** s1 = data[1];
+
+  if (gts_bbox_overlaps_triangle (bb, t)) {
+    if (*s1 == NULL)
+      *s1 = gts_surface_new (gts_surface_class (),
+			     gts_face_class (),
+			     gts_edge_class (),
+			     gts_vertex_class ());
+    gts_surface_add_face (*s1, GTS_FACE (t));
+  }
+}
+
+/**
+ * gfs_cell_is_cut:
+ * @cell: a #FttCell.
+ * @s: a #GfsSurface.
+ * @flatten: if set to %TRUE, @cell is flattened in the z direction.
+ *
+ * Returns: a new #GfsSurface containing a subset of @s which may
+ * intersect @cell or %NULL if @s does not intersect @cell.
+ */
+GfsSurface * gfs_cell_is_cut (FttCell * cell, GfsSurface * s, gboolean flatten)
+{
+  g_return_val_if_fail (cell != NULL, NULL);
+  g_return_val_if_fail (s != NULL, NULL);
+
+  if (s->s) {
+    GtsSurface * s1 = NULL;
+    gpointer data[2];
+    GtsBBox bb;
+    
+    ftt_cell_bbox (cell, &bb);
+    if (flatten)
+      bb.z1 = bb.z2 = 0.;
+    data[0] = &bb;
+    data[1] = &s1;
+    gts_surface_foreach_face (s->s, (GtsFunc) face_overlaps_box, data);
+    if (s1 == NULL)
+      return NULL;
+    GfsSurface * s2 = GFS_SURFACE (gts_object_new (gfs_surface_class ()));
+    s2->s = s1;
+    return s2;
+  }
+
+  g_assert_not_implemented ();
+}
+
+static void cell_traverse_cut (FttCell * cell,
+			       GfsSurface * s,
+			       FttTraverseType order,
+			       FttTraverseFlags flags,
+			       FttCellTraverseCutFunc func,
+			       gpointer data,
+			       gboolean flatten)
+{
+  GfsSurface * s1 = gfs_cell_is_cut (cell, s, flatten && FTT_CELL_IS_LEAF (cell));
+
+  if (s1 == NULL)
+    return;
+  if (order == FTT_PRE_ORDER &&
+      (flags == FTT_TRAVERSE_ALL ||
+       ((flags & FTT_TRAVERSE_LEAFS) != 0 && FTT_CELL_IS_LEAF (cell)) ||
+       ((flags & FTT_TRAVERSE_NON_LEAFS) != 0 && !FTT_CELL_IS_LEAF (cell))))
+    (* func) (cell, s1, data);
+  if (!FTT_CELL_IS_LEAF (cell)) {
+    struct _FttOct * children = cell->children;
+    guint n;
+
+    for (n = 0; n < FTT_CELLS; n++) {
+      FttCell * c = &(children->cell[n]);
+
+      if (!FTT_CELL_IS_DESTROYED (c))
+	cell_traverse_cut (c, s1, order, flags, func, data, flatten);
+    }
+  }
+  if (order == FTT_POST_ORDER &&
+      (flags == FTT_TRAVERSE_ALL ||
+       ((flags & FTT_TRAVERSE_LEAFS) != 0 && FTT_CELL_IS_LEAF (cell)) ||
+       ((flags & FTT_TRAVERSE_NON_LEAFS) != 0 && !FTT_CELL_IS_LEAF (cell))))
+    (* func) (cell, s1, data);
+  gts_object_destroy (GTS_OBJECT (s1));
+}
+
+/**
+ * gfs_cell_traverse_cut:
+ * @root: the root #FttCell of the tree to traverse.
+ * @s: a #GfsSurface.
+ * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
+ * %FTT_POST_ORDER. 
+ * @flags: which types of children are to be visited.
+ * @func: the function to call for each visited #FttCell.
+ * @data: user data to pass to @func.
+ * 
+ * Traverses a cell tree starting at the given root #FttCell. Calls
+ * the given function for each cell cut by @s.
+ */
+void gfs_cell_traverse_cut (FttCell * root,
+			    GfsSurface * s,
+			    FttTraverseType order,
+			    FttTraverseFlags flags,
+			    FttCellTraverseCutFunc func,
+			    gpointer data)
+{
+  g_return_if_fail (root != NULL);
+  g_return_if_fail (s != NULL);
+  g_return_if_fail (func != NULL);
+
+  cell_traverse_cut (root, s, order, flags, func, data, FALSE);
+}
+
+/**
+ * gfs_cell_traverse_cut_2D:
+ * @root: the root #FttCell of the tree to traverse.
+ * @s: a #GfsSurface.
+ * @order: the order in which the cells are visited - %FTT_PRE_ORDER,
+ * %FTT_POST_ORDER. 
+ * @flags: which types of children are to be visited.
+ * @func: the function to call for each visited #FttCell.
+ * @data: user data to pass to @func.
+ * 
+ * Traverses a cell tree starting at the given root #FttCell. Calls
+ * the given function for each cell cut by @s.
+ *
+ * The cells are "flattened" in the z-direction.
+ */
+void gfs_cell_traverse_cut_2D (FttCell * root,
+			       GfsSurface * s,
+			       FttTraverseType order,
+			       FttTraverseFlags flags,
+			       FttCellTraverseCutFunc func,
+			       gpointer data)
+{
+  g_return_if_fail (root != NULL);
+  g_return_if_fail (s != NULL);
+  g_return_if_fail (func != NULL);
+
+  cell_traverse_cut (root, s, order, flags, func, data, TRUE);
+}
diff --git a/src/surface.h b/src/surface.h
new file mode 100644
index 0000000..13ab673
--- /dev/null
+++ b/src/surface.h
@@ -0,0 +1,73 @@
+/* Gerris - The GNU Flow Solver
+ * Copyright (C) 2001 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 __SURFACE_H__
+#define __SURFACE_H__
+
+#include <gts.h>
+#include "ftt.h"
+
+/* GfsSurface: Header */
+
+typedef struct _GfsSurface         GfsSurface;
+
+struct _GfsSurface {
+  /*< private >*/
+  GtsObject parent;
+
+  /*< public >*/
+  GtsSurface * s;
+};
+
+#define GFS_SURFACE(obj)            GTS_OBJECT_CAST (obj,\
+					         GfsSurface,\
+					         gfs_surface_class ())
+#define GFS_IS_SURFACE(obj)         (gts_object_is_from_class (obj,\
+						 gfs_surface_class ()))
+
+GtsObjectClass *   gfs_surface_class          (void);
+void               gfs_surface_read           (GfsSurface * s, 
+					       GtsFile * fp);
+void               gfs_surface_write          (GfsSurface * s,
+					       gpointer sim,
+					       FILE * fp);
+GfsSurface *       gfs_cell_is_cut            (FttCell * cell,
+					       GfsSurface * s,
+					       gboolean flatten);
+typedef void       (* FttCellTraverseCutFunc) (FttCell * cell,
+					       GfsSurface * s,
+					       gpointer data);
+void               gfs_cell_traverse_cut       (FttCell * root,
+						GfsSurface * s,
+						FttTraverseType order,
+						FttTraverseFlags flags,
+						FttCellTraverseCutFunc func,
+						gpointer data);
+void               gfs_cell_traverse_cut_2D    (FttCell * root,
+						GfsSurface * s,
+						FttTraverseType order,
+						FttTraverseFlags flags,
+						FttCellTraverseCutFunc func,
+						gpointer data);
+
+#ifdef __cplusplus
+}
+#endif /* __cplusplus */
+
+#endif /* __SURFACE_H__ */
diff --git a/src/utils.c b/src/utils.c
index 6470d72..c7b7af9 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -1008,7 +1008,7 @@ GtsObjectClass * gfs_object_class_from_name (const gchar * name)
     return klass;
   /* for backward parameter file compatibility */
   if (!strcmp (name, "GtsSurfaceFile"))
-    return GTS_OBJECT_CLASS (gfs_surface_class ());
+    return GTS_OBJECT_CLASS (gfs_solid_class ());
   gchar * ename = g_strconcat ("Gfs", name, NULL);
   klass = gts_object_class_from_name (ename);
   g_free (ename);
diff --git a/test/boundaries/boundaries.gfs b/test/boundaries/boundaries.gfs
index 67e2162..cc3be9d 100644
--- a/test/boundaries/boundaries.gfs
+++ b/test/boundaries/boundaries.gfs
@@ -34,7 +34,7 @@
   AdvectionParams { scheme = none }
   ApproxProjectionParams { tolerance = 1e-6 }
   Refine LEVEL
-  Surface boundaries.gts
+  Solid boundaries.gts
   Init {} { U = 1 }
   OutputSimulation { start = end } sim-LEVEL {
     variables = U,V,P
diff --git a/test/channel/channel.gfs b/test/channel/channel.gfs
index ceafa11..c60e8c7 100644
--- a/test/channel/channel.gfs
+++ b/test/channel/channel.gfs
@@ -34,7 +34,7 @@
   ProjectionParams { tolerance = 1e-6 }
   ApproxProjectionParams { tolerance = 1e-6 }
   Refine LEVEL
-  Surface channel.gts
+  Solid channel.gts
   Init {} { U = 1 }
   OutputSimulation { start = end } sim-LEVEL {
     variables = U,V,P
diff --git a/test/circle/circle.gfs b/test/circle/circle.gfs
index 548ba97..cbcf268 100644
--- a/test/circle/circle.gfs
+++ b/test/circle/circle.gfs
@@ -54,7 +54,7 @@
 1 0 GfsPoisson GfsBox GfsGEdge {} {
   Time { iend = 10 }
   Refine LEVEL
-  Surface shape.gts
+  Solid shape.gts
   ApproxProjectionParams { nrelax = 4 tolerance = 1e-30 erelax = 2 }
   Init {} {
     Div = {
diff --git a/test/circle/refined/refined.gfs b/test/circle/refined/refined.gfs
index 8b627c0..4f1435e 100644
--- a/test/circle/refined/refined.gfs
+++ b/test/circle/refined/refined.gfs
@@ -50,7 +50,7 @@
   Time { iend = 10 }
   Refine LEVEL
   RefineSolid (LEVEL + 2)
-  Surface shape.gts
+  Solid shape.gts
   ApproxProjectionParams { nrelax = 4 tolerance = 1e-30 }
   Init {} {
     Div = {
diff --git a/test/circle/star/star.gfs b/test/circle/star/star.gfs
index 51baa81..d4b82e3 100644
--- a/test/circle/star/star.gfs
+++ b/test/circle/star/star.gfs
@@ -56,7 +56,7 @@
 1 0 GfsPoisson GfsBox GfsGEdge {} {
   Time { iend = 10 }
   Refine LEVEL
-  Surface shape.gts
+  Solid shape.gts
   ApproxProjectionParams { nrelax = 4 tolerance = 1e-30 erelax = 2 }
   Init {} {
     Div = {
diff --git a/test/couette/couette.gfs b/test/couette/couette.gfs
index ef45735..7f46101 100644
--- a/test/couette/couette.gfs
+++ b/test/couette/couette.gfs
@@ -40,8 +40,8 @@
 1 0 GfsSimulation GfsBox GfsGEdge {} {
   Time { iend = 100 dtmax = 1e-2 }
   Refine 6
-  Surface inner.gts { flip = 1 scale = 1.9999 }
-  Surface inner.gts
+  Solid inner.gts { flip = 1 scale = 1.9999 }
+  Solid inner.gts
   ApproxProjectionParams { tolerance = 1e-6 }
   AdvectionParams { scheme = none }
   SourceViscosity {} {
diff --git a/test/dumbell/dumbell.gfs b/test/dumbell/dumbell.gfs
index f1ccf40..155a9c3 100644
--- a/test/dumbell/dumbell.gfs
+++ b/test/dumbell/dumbell.gfs
@@ -13,7 +13,7 @@
 1 0 GfsPoisson GfsBox GfsGEdge {} {
   Refine 3
   ApproxProjectionParams { nitermax = 1000 minlevel = 1 tolerance = 1e-30 }
-  Surface dumbell.gts
+  Solid dumbell.gts
   Init {} { Div = y }
   OutputProjectionStats { istep = 1 } stdout
 }
diff --git a/test/geo/geo.gfs b/test/geo/geo.gfs
index 0bc5749..c65e5da 100644
--- a/test/geo/geo.gfs
+++ b/test/geo/geo.gfs
@@ -98,9 +98,9 @@
     unbiased = 1
     v = E
   }
-  OutputSimulation { istart = 100 iend = 500 istep = 100 } stdout {}
+  OutputSimulation { istart = 100 iend = 500 istep = 100 } stdout
   EventScript { istart = 100 iend = 500 istep = 100 } { echo "Save error-$GfsIter.eps { format = EPS }"}
-  OutputSimulation { istart = 1500 } stdout {}
+  OutputSimulation { istart = 1500 } stdout
   EventScript { istart = 1500 } { echo "Save error-$GfsIter.eps { format = EPS }"}
   EventScript { start = end } {
     cat <<EOF | gnuplot
diff --git a/test/merging/merging.gfs b/test/merging/merging.gfs
index 40dd168..a09b598 100644
--- a/test/merging/merging.gfs
+++ b/test/merging/merging.gfs
@@ -121,7 +121,7 @@
            vortex (-0.045, -0.0779422863406, 50.);
   }
   AdaptVorticity { istep = 1 } { maxlevel = LEVEL cmax = 4e-3 }
-  OutputSimulation { start = 0.05 } stdout {}
+  OutputSimulation { start = 0.05 } stdout
   EventScript { start = 0.05 } {
     echo Clear
     cat levels.gfv
@@ -130,7 +130,7 @@
     cat vorticity.gfv
     echo Save tv_0_05.eps { format = EPS line_width = 0.1 }
   }
-  OutputSimulation { start = 0.15 } stdout {}
+  OutputSimulation { start = 0.15 } stdout
   EventScript { start = 0.15 } {
     echo Clear
     cat levels.gfv
@@ -139,7 +139,7 @@
     cat vorticity.gfv
     echo Save tv_0_15.eps { format = EPS line_width = 0.1 }
   }
-  OutputSimulation { start = 0.25 } stdout {}
+  OutputSimulation { start = 0.25 } stdout
   EventScript { start = 0.25 } {
     echo Clear
     cat levels.gfv
@@ -148,6 +148,6 @@
     cat vorticity.gfv
     echo Save tv_0_25.eps { format = EPS line_width = 0.1 }
   }
-  OutputSimulation { start = 0.25 } SIM {}
+  OutputSimulation { start = 0.25 } SIM
 }
 GfsBox {}
diff --git a/test/nz/nz.gfs b/test/nz/nz.gfs
index 3e710d3..ce10668 100644
--- a/test/nz/nz.gfs
+++ b/test/nz/nz.gfs
@@ -45,7 +45,7 @@
     AdvectionParams { scheme = none }
     Init {} { P = 1e-2*x }
     # Bathymetry
-    Surface bath.gts
+    Solid bath.gts
     # Refine the coastline to 7 levels
     RefineSurface 7 bath.gts { twod = 1 }
     # Non-dimensional gravity
diff --git a/test/plate/plate.gfs b/test/plate/plate.gfs
index bb25f33..ed9d5ea 100644
--- a/test/plate/plate.gfs
+++ b/test/plate/plate.gfs
@@ -16,7 +16,7 @@
   Time { iend = 30 dtmax = 1e-2 }
   Refine 5
   RefineSolid 6
-  Surface square.gts
+  Solid square.gts
   AdvectionParams { scheme = none }
   Init {} { U = 1 }
   OutputScalarNorm { start = end } stdout { v = Velocity } 
diff --git a/test/waves/adaptive/adaptive.gfs b/test/waves/adaptive/adaptive.gfs
index ecb60c3..d3b29e1 100644
--- a/test/waves/adaptive/adaptive.gfs
+++ b/test/waves/adaptive/adaptive.gfs
@@ -81,7 +81,7 @@
   }
   Refine LEVEL
   AdaptVorticity { istep = 1 } { cmax = 5e-2 maxlevel = LEVEL }
-  Surface basin.gts
+  Solid basin.gts
   AdvectionParams { scheme = none }
   ApproxProjectionParams { tolerance = 1e-9 weighted = 0 }
   SourceCoriolis {} 1.
@@ -617,7 +617,7 @@
     s = pwave(cx, cy, 30)
     unbiased = 1
   }
-  OutputSimulation { start = end } sim-LEVEL {}
+  OutputSimulation { start = end } sim-LEVEL
 }
 GfsBox {
   front = Boundary
diff --git a/test/waves/waves.gfs b/test/waves/waves.gfs
index f198713..24c6128 100644
--- a/test/waves/waves.gfs
+++ b/test/waves/waves.gfs
@@ -92,7 +92,7 @@
     H = 1
   }
   Refine LEVEL
-  Surface basin.gts
+  Solid basin.gts
   AdvectionParams { scheme = none }
   ApproxProjectionParams { tolerance = 1e-9 weighted = 0 }
   SourceCoriolis {} 1.
@@ -628,7 +628,7 @@
     s = pwave(cx, cy, 30)
     unbiased = 1
   }
-  OutputSimulation { start = end } sim-LEVEL {}
+  OutputSimulation { start = end } sim-LEVEL
 }
 GfsBox {
   front = Boundary
diff --git a/tools/gfs2gfs b/tools/gfs2gfs
index 849f96a..5e49838 100755
--- a/tools/gfs2gfs
+++ b/tools/gfs2gfs
@@ -30,5 +30,6 @@ while test $# -gt 0; do
   shift
 done
 
-sed 's/^ *GtsSurface/GfsSurface {}/g' | \
-sed 's/GtsSurfaceFile/GfsSurface/g'
+sed 's/^ *GtsSurface/GfsSolid {}/g' | \
+sed 's/GtsSurfaceFile/GfsSolid/g' | \
+sed 's/surface =/solid =/g'
diff --git a/tools/gfs2oogl.c b/tools/gfs2oogl.c
index 69f6ab5..7a8b64c 100644
--- a/tools/gfs2oogl.c
+++ b/tools/gfs2oogl.c
@@ -1154,21 +1154,24 @@ int main (int argc, char * argv[])
       printf ("})\n");
     }
     else if (draw_surface) {
-      GtsSurface * s = surface;
-
-      if (s == NULL)
-	s = gfs_simulation_get_surface (simulation);
-
-      if (s) {
-	if (refine)
-	  gts_surface_refine (s, 
-			      (GtsKeyFunc) local_size_ratio, domain,
-			      NULL, NULL,
-			      (GtsStopFunc) stop, NULL);
-	gfs_draw_surface (domain, s, 
-			  var, stats.min, stats.max,
-			  stdout);
+      GSList * l = gfs_simulation_get_solids (simulation), * i = l;
+	
+      while (i) {
+	GfsSurface * s = i->data;
+	  
+	if (s->s) {
+	  if (refine)
+	    gts_surface_refine (s->s, 
+				(GtsKeyFunc) local_size_ratio, domain,
+				NULL, NULL,
+				(GtsStopFunc) stop, NULL);
+	  gfs_draw_surface (domain, s->s, 
+			    var, stats.min, stats.max,
+			    stdout);
+	}
+	i = i->next;
       }
+      g_slist_free (l);
     }
     else if (merged) {
       gfs_set_merged (domain);

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list