diff --git a/Plugin/Annotate.cpp b/Plugin/Annotate.cpp
index ddd078506cfd7bb65f309db5beda1030d9f37857..ae05e44c16abc161152a3125678963b8c9189d23 100644
--- a/Plugin/Annotate.cpp
+++ b/Plugin/Annotate.cpp
@@ -1,4 +1,4 @@
-// $Id: Annotate.cpp,v 1.4 2005-01-01 19:35:37 geuzaine Exp $
+// $Id: Annotate.cpp,v 1.5 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -73,14 +73,13 @@ void GMSH_AnnotatePlugin::getInfos(char *author, char *copyright,
   strcpy(copyright, "DGR (www.multiphysics.com)");
   strcpy(help_text,
          "Plugin(Annotate) adds a text string of size\n"
-	 "`FontSize' in the view `iView'. If `3D' is\n"
-	 "equal to 1, the plugin inserts the string\n"
-	 "in model coordinates at the position (`X',`Y',`Z').\n"
-	 "If `3D' is equal to 0, the plugin inserts the\n"
-	 "string in screen coordinates at the position\n"
-	 "(`X',`Y'), and aligns it according to `Align'.\n"
-	 "If `iView' < 0, the plugin is run on the current\n"
-	 "view.\n"
+	 "`FontSize' in the view `iView'. If `3D' is equal\n"
+	 "to 1, the plugin inserts the string in model\n"
+	 "coordinates at the position (`X',`Y',`Z'). If\n"
+	 "`3D' is equal to 0, the plugin inserts the string\n"
+	 "in screen coordinates at the position (`X',`Y'),\n"
+	 "and aligns it according to `Align'. If `iView'\n"
+	 "< 0, the plugin is run on the current view.\n"
 	 "\n"
 	 "Plugin(Annotate) is executed in-place.\n");
 }
diff --git a/Plugin/CutMap.cpp b/Plugin/CutMap.cpp
index a3a582bc6bc5bd891e12432c2a3713124b8ce9f7..e2ae355cacd9b4acd28541912f8743eb30e163d7 100644
--- a/Plugin/CutMap.cpp
+++ b/Plugin/CutMap.cpp
@@ -1,4 +1,4 @@
-// $Id: CutMap.cpp,v 1.46 2005-01-01 19:35:37 geuzaine Exp $
+// $Id: CutMap.cpp,v 1.47 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -29,9 +29,10 @@ StringXNumber CutMapOptions_Number[] = {
   {GMSH_FULLRC, "A", GMSH_CutMapPlugin::callbackA, 1.},
   {GMSH_FULLRC, "dTimeStep", NULL, -1.},
   {GMSH_FULLRC, "dView", NULL, -1.},
-  {GMSH_FULLRC, "iView", NULL, -1.},
-  {GMSH_FULLRC, "recurLevel", NULL, 4},
-  {GMSH_FULLRC, "targetError", NULL, 0}
+  {GMSH_FULLRC, "extractVolume", GMSH_CutMapPlugin::callbackVol, 0.},
+  {GMSH_FULLRC, "recurLevel", GMSH_CutMapPlugin::callbackRecur, 4},
+  {GMSH_FULLRC, "targetError", GMSH_CutMapPlugin::callbackTarget, 0},
+  {GMSH_FULLRC, "iView", NULL, -1.}
 };
 
 extern "C"
@@ -51,7 +52,7 @@ double GMSH_CutMapPlugin::callbackA(int num, int action, double value)
 {
   double min = 0., max = 1.;
   if(action > 0){
-    int iview = (int)CutMapOptions_Number[3].def;
+    int iview = (int)CutMapOptions_Number[6].def;
     if(iview < 0) iview = num;
     Post_View **vv = (Post_View **)List_Pointer_Test(CTX.post.list, iview);
     if(vv){
@@ -68,6 +69,39 @@ double GMSH_CutMapPlugin::callbackA(int num, int action, double value)
   return 0.;
 }
 
+double GMSH_CutMapPlugin::callbackVol(int num, int action, double value)
+{
+  switch(action){ // configure the input field
+  case 1: return 1.;
+  case 2: return -1.;
+  case 3: return 1.;
+  default: break;
+  }
+  return 0.;
+}
+
+double GMSH_CutMapPlugin::callbackRecur(int num, int action, double value)
+{
+  switch(action){ // configure the input field
+  case 1: return 1.;
+  case 2: return 0.;
+  case 3: return 10.;
+  default: break;
+  }
+  return 0.;
+}
+
+double GMSH_CutMapPlugin::callbackTarget(int num, int action, double value)
+{
+  switch(action){ // configure the input field
+  case 1: return 0.01;
+  case 2: return 0.;
+  case 3: return 1.;
+  default: break;
+  }
+  return 0.;
+}
+
 void GMSH_CutMapPlugin::getName(char *name) const
 {
   strcpy(name, "Cut Map");
@@ -79,14 +113,18 @@ void GMSH_CutMapPlugin::getInfos(char *author, char *copyright,
   strcpy(author, "J.-F. Remacle (remacle@scorec.rpi.edu)");
   strcpy(copyright, "DGR (www.multiphysics.com)");
   strcpy(help_text,
-         "Plugin(CutMap) extracts the isovalue surface of\n"
-         "value `A' from the view `iView' and draws the\n"
+         "Plugin(CutMap) extracts the isosurface of value\n"
+         "`A' from the view `iView' and draws the\n"
 	 "`dTimeStep'-th value of the view `dView' on this\n"
 	 "isovalue surface. If `iView' < 0, the plugin is\n"
 	 "run on the current view. If `dTimeStep' < 0, the\n"
 	 "plugin uses, for each time step in `iView', the\n"
 	 "corresponding time step in `dView'. If `dView'\n"
 	 "< 0, the plugin uses `iView' as the value source.\n"
+	 "If `extractVolume' is nonzero, the plugin\n" 
+	 "extracts the isovolume with values smaller (if\n"
+	 "`extractVolume' > 0) or greater (if `extractVolume'\n"
+	 "< 0) than the isosurface `A'.\n"
 	 "\n"
 	 "Plugin(CutMap) creates as many views as there\n"
 	 "are time steps in `iView'.\n");
@@ -117,10 +155,11 @@ double GMSH_CutMapPlugin::levelset(double x, double y, double z, double val) con
 
 Post_View *GMSH_CutMapPlugin::execute(Post_View * v)
 {
-  int iView = (int)CutMapOptions_Number[3].def;
+  int iView = (int)CutMapOptions_Number[6].def;
   _valueIndependent = 0;
   _valueView = (int)CutMapOptions_Number[2].def;
   _valueTimeStep = (int)CutMapOptions_Number[1].def;
+  _extractVolume = (int)CutMapOptions_Number[3].def;
   _recurLevel = (int)CutMapOptions_Number[4].def;
   _targetError = CutMapOptions_Number[5].def;
   _orientation = GMSH_LevelsetPlugin::MAP;
diff --git a/Plugin/CutMap.h b/Plugin/CutMap.h
index 0ea6103b0186fce381ea1d5e4e02588d2dab0df7..5a53fb7d898e42c7c874d0233600dd23b1d70321 100644
--- a/Plugin/CutMap.h
+++ b/Plugin/CutMap.h
@@ -39,6 +39,9 @@ class GMSH_CutMapPlugin : public GMSH_LevelsetPlugin
   StringXNumber* getOption (int iopt);  
   Post_View *execute (Post_View *);
   static double callbackA(int, int, double);
+  static double callbackVol(int, int, double);
+  static double callbackRecur(int, int, double);
+  static double callbackTarget(int, int, double);
 };
 
 #endif
diff --git a/Plugin/CutPlane.cpp b/Plugin/CutPlane.cpp
index 244d4ef8030eedf36a584a709b17c2030fd06f1b..30e7b20ce70d8fe35c1d4ae62ba6c3f1becbc56d 100644
--- a/Plugin/CutPlane.cpp
+++ b/Plugin/CutPlane.cpp
@@ -1,4 +1,4 @@
-// $Id: CutPlane.cpp,v 1.44 2005-01-01 19:35:38 geuzaine Exp $
+// $Id: CutPlane.cpp,v 1.45 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -37,9 +37,10 @@ StringXNumber CutPlaneOptions_Number[] = {
   {GMSH_FULLRC, "B", GMSH_CutPlanePlugin::callbackB, 0.},
   {GMSH_FULLRC, "C", GMSH_CutPlanePlugin::callbackC, 0.},
   {GMSH_FULLRC, "D", GMSH_CutPlanePlugin::callbackD, -0.01},
-  {GMSH_FULLRC, "iView", NULL, -1.},
-  {GMSH_FULLRC, "recurLevel", NULL, 4},
-  {GMSH_FULLRC, "targetError", NULL, 0}
+  {GMSH_FULLRC, "extractVolume", GMSH_CutPlanePlugin::callbackVol, 0},
+  {GMSH_FULLRC, "recurLevel", GMSH_CutPlanePlugin::callbackRecur, 4},
+  {GMSH_FULLRC, "targetError", GMSH_CutPlanePlugin::callbackTarget, 0.},
+  {GMSH_FULLRC, "iView", NULL, -1.}
 };
 
 extern "C"
@@ -58,7 +59,7 @@ GMSH_CutPlanePlugin::GMSH_CutPlanePlugin()
 void GMSH_CutPlanePlugin::draw()
 {
 #if defined(HAVE_FLTK)
-  int num = (int)CutPlaneOptions_Number[4].def;
+  int num = (int)CutPlaneOptions_Number[7].def;
   if(num < 0) num = iview;
   Post_View **vv = (Post_View **)List_Pointer_Test(CTX.post.list, num);
   if(!vv) return;
@@ -113,6 +114,24 @@ double GMSH_CutPlanePlugin::callbackD(int num, int action, double value)
 		  CTX.lc/200., -CTX.lc, CTX.lc);
 }
 
+double GMSH_CutPlanePlugin::callbackVol(int num, int action, double value)
+{
+  return callback(num, action, value, &CutPlaneOptions_Number[4].def,
+		  1., -1, 1);
+}
+
+double GMSH_CutPlanePlugin::callbackRecur(int num, int action, double value)
+{
+  return callback(num, action, value, &CutPlaneOptions_Number[5].def,
+		  1, 0, 10);
+}
+
+double GMSH_CutPlanePlugin::callbackTarget(int num, int action, double value)
+{
+  return callback(num, action, value, &CutPlaneOptions_Number[6].def,
+		  0.01, 0., 1.);
+}
+
 void GMSH_CutPlanePlugin::getName(char *name) const
 {
   strcpy(name, "Cut Plane");
@@ -126,8 +145,10 @@ void GMSH_CutPlanePlugin::getInfos(char *author, char *copyright,
   strcpy(help_text,
          "Plugin(CutPlane) cuts the view `iView' with\n"
 	 "the plane `A'*X + `B'*Y + `C'*Z + `D' = 0. If\n"
-	 "`iView' < 0, the plugin is run on the current\n"
-	 "view.\n"
+	 "`extractVolume' is nonzero, the plugin extracts\n"
+	 "the elements on one side of the plane (depending\n"
+	 "on the sign of `extractVolume'). If `iView' < 0,\n"
+	 "the plugin is run on the current view.\n"
 	 "\n"
 	 "Plugin(CutPlane) creates one new view.\n");
 }
@@ -156,7 +177,7 @@ double GMSH_CutPlanePlugin::levelset(double x, double y, double z, double val) c
 
 Post_View *GMSH_CutPlanePlugin::execute(Post_View * v)
 {
-  int iView = (int)CutPlaneOptions_Number[4].def;
+  int iView = (int)CutPlaneOptions_Number[7].def;
   _ref[0] = CutPlaneOptions_Number[0].def;
   _ref[1] = CutPlaneOptions_Number[1].def;
   _ref[2] = CutPlaneOptions_Number[2].def;
@@ -164,6 +185,7 @@ Post_View *GMSH_CutPlanePlugin::execute(Post_View * v)
   _valueView = -1;
   _valueTimeStep = -1;
   _orientation = GMSH_LevelsetPlugin::PLANE;
+  _extractVolume = (int)CutPlaneOptions_Number[4].def;
   _recurLevel = (int)CutPlaneOptions_Number[5].def;
   _targetError = CutPlaneOptions_Number[6].def;
   
diff --git a/Plugin/CutPlane.h b/Plugin/CutPlane.h
index 3546f6a4a8ee76aa0f3a1efd0d3daec8e85f5c7d..a8b0ed97db7b328ed20a03236f0babd02b90f5b0 100644
--- a/Plugin/CutPlane.h
+++ b/Plugin/CutPlane.h
@@ -46,6 +46,9 @@ public:
   static double callbackB(int, int, double);
   static double callbackC(int, int, double);
   static double callbackD(int, int, double);
+  static double callbackVol(int, int, double);
+  static double callbackRecur(int, int, double);
+  static double callbackTarget(int, int, double);
   static void draw();
 };
 
diff --git a/Plugin/CutSphere.cpp b/Plugin/CutSphere.cpp
index c6b8f4b433a447ce07b61bdfc88fc7e37fad282c..5683a47ec97219c6f96a95fb61beacee4b3c1cfd 100644
--- a/Plugin/CutSphere.cpp
+++ b/Plugin/CutSphere.cpp
@@ -1,4 +1,4 @@
-// $Id: CutSphere.cpp,v 1.41 2005-01-01 19:35:38 geuzaine Exp $
+// $Id: CutSphere.cpp,v 1.42 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -36,8 +36,9 @@ StringXNumber CutSphereOptions_Number[] = {
   {GMSH_FULLRC, "Yc", GMSH_CutSpherePlugin::callbackY, 0.},
   {GMSH_FULLRC, "Zc", GMSH_CutSpherePlugin::callbackZ, 0.},
   {GMSH_FULLRC, "R", GMSH_CutSpherePlugin::callbackR, 0.25},
-  {GMSH_FULLRC, "iView", NULL, -1.},
-  {GMSH_FULLRC, "recurLevel", NULL, 4}
+  {GMSH_FULLRC, "extractVolume", GMSH_CutSpherePlugin::callbackVol, 0.},
+  {GMSH_FULLRC, "recurLevel", GMSH_CutSpherePlugin::callbackRecur, 4},
+  {GMSH_FULLRC, "iView", NULL, -1.}
 };
 
 extern "C"
@@ -117,6 +118,18 @@ double GMSH_CutSpherePlugin::callbackR(int num, int action, double value)
 		  CTX.lc/200., 0., 2 * CTX.lc);
 }
 
+double GMSH_CutSpherePlugin::callbackVol(int num, int action, double value)
+{
+  return callback(num, action, value, &CutSphereOptions_Number[4].def,
+		  1., -1., 1.);
+}
+
+double GMSH_CutSpherePlugin::callbackRecur(int num, int action, double value)
+{
+  return callback(num, action, value, &CutSphereOptions_Number[5].def,
+		  1, 0, 10);
+}
+
 void GMSH_CutSpherePlugin::getName(char *name) const
 {
   strcpy(name, "Cut Sphere");
@@ -130,8 +143,10 @@ void GMSH_CutSpherePlugin::getInfos(char *author, char *copyright,
   strcpy(help_text,
          "Plugin(CutSphere) cuts the view `iView' with the\n"
 	 "sphere (X-`Xc')^2 + (Y-`Yc')^2 + (Z-`Zc')^2 = `R'^2.\n"
-	 "If `iView' < 0, the plugin is run on the current\n"
-	 "view.\n"
+	 "If `extractVolume' is nonzero, the plugin extracts\n"
+	 "the elements inside (if `extractVolume' < 0) or\n"
+	 "outside (if `extractVolume' > 0) the sphere. If\n"
+	 "`iView' < 0, the plugin is run on the current view.\n"
 	 "\n"
 	 "Plugin(CutSphere) creates one new view.\n");
 }
@@ -163,10 +178,11 @@ double GMSH_CutSpherePlugin::levelset(double x, double y, double z,
 
 Post_View *GMSH_CutSpherePlugin::execute(Post_View * v)
 {
-  int iView = (int)CutSphereOptions_Number[4].def;
+  int iView = (int)CutSphereOptions_Number[6].def;
   _ref[0] = CutSphereOptions_Number[0].def;
   _ref[1] = CutSphereOptions_Number[1].def;
   _ref[2] = CutSphereOptions_Number[2].def;
+  _extractVolume = (int)CutSphereOptions_Number[4].def;
   _recurLevel = (int)CutSphereOptions_Number[5].def;
 
   _valueIndependent = 1;
diff --git a/Plugin/CutSphere.h b/Plugin/CutSphere.h
index 3d147e1289dbea50ff0aaf8184a52a49ef380741..ca7e808721e61554e5f433ab7fcfbaeae4c4d31c 100644
--- a/Plugin/CutSphere.h
+++ b/Plugin/CutSphere.h
@@ -45,6 +45,8 @@ public:
   static double callbackY(int, int, double);
   static double callbackZ(int, int, double);
   static double callbackR(int, int, double);
+  static double callbackVol(int, int, double);
+  static double callbackRecur(int, int, double);
   static void draw();
 };
 
diff --git a/Plugin/DecomposeInSimplex.cpp b/Plugin/DecomposeInSimplex.cpp
index bf564e5f4866c9ebc8bc7cc1451e1c3a02e371a7..31dd8b2f66703b1954f30e5f22e46cef934ecb94 100644
--- a/Plugin/DecomposeInSimplex.cpp
+++ b/Plugin/DecomposeInSimplex.cpp
@@ -1,4 +1,4 @@
-// $Id: DecomposeInSimplex.cpp,v 1.16 2005-01-01 19:35:38 geuzaine Exp $
+// $Id: DecomposeInSimplex.cpp,v 1.17 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -26,6 +26,7 @@
 #include "Views.h"
 #include "Context.h"
 #include "Malloc.h"
+#include "Numeric.h"
 
 extern Context_T CTX;
 
@@ -200,11 +201,32 @@ void DecomposeInSimplex::reorder(int map[4], int n,
     zn[i] = z[map[i]];
   }
 
+  int map2[4] = {map[0], map[1], map[2], map[3]};
+#if 0
+  // make sure to create tets with positive volume?
+  if(n == 4){ // tets
+    double mat[3][3];
+    mat[0][0] = xn[1] - xn[0]; mat[0][1] = xn[2] - xn[0]; mat[0][2] = xn[3] - xn[0];
+    mat[1][0] = yn[1] - yn[0]; mat[1][1] = yn[2] - yn[0]; mat[1][2] = yn[3] - yn[0];
+    mat[2][0] = zn[1] - zn[0]; mat[2][1] = zn[2] - zn[0]; mat[2][2] = zn[3] - zn[0];
+    if(det3x3(mat) < 0.){
+      map2[0] = map[1];
+      map2[1] = map[0];
+      xn[0] = x[map2[0]];
+      yn[0] = y[map2[0]];
+      zn[0] = z[map2[0]];
+      xn[1] = x[map2[1]];
+      yn[1] = y[map2[1]];
+      zn[1] = z[map2[1]];
+    }
+  }
+#endif
+
   for(int ts = 0; ts < _numTimeSteps; ts++)
     for(int i = 0; i < n; i++) {
       for(int j = 0; j < _numComponents; j++)
 	valn[ts*n*_numComponents + i*_numComponents + j] = 
-	  val[ts*_numNodes*_numComponents + map[i]*_numComponents + j];
+	  val[ts*_numNodes*_numComponents + map2[i]*_numComponents + j];
   }
 }
 
diff --git a/Plugin/DisplacementRaise.cpp b/Plugin/DisplacementRaise.cpp
index 630f7843bebd3afe2c52bcc219f1834750fbd162..1f9bb3fe1979897188b2455be50d2c3c261b90d7 100644
--- a/Plugin/DisplacementRaise.cpp
+++ b/Plugin/DisplacementRaise.cpp
@@ -1,4 +1,4 @@
-// $Id: DisplacementRaise.cpp,v 1.19 2005-01-01 19:35:39 geuzaine Exp $
+// $Id: DisplacementRaise.cpp,v 1.20 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -64,11 +64,10 @@ void GMSH_DisplacementRaisePlugin::getInfos(char *author, char *copyright,
          "coordinates of the elements in the view `iView'\n"
          "using the vectorial values (the displacements)\n"
          "stored in the `dTimeStep'-th time step of the\n"
-	 "view `dView'. If `iView' < 0, the plugin is\n"
-	 "run on the current view. If `dView' < 0, the\n"
-	 "plugin looks for the displacements in the\n"
-	 "view located just after `iView' in the view\n"
-	 "list.\n"
+	 "view `dView'. If `iView' < 0, the plugin is run\n"
+	 "on the current view. If `dView' < 0, the plugin\n"
+	 "looks for the displacements in the view located\n"
+	 "just after `iView' in the view list.\n"
 	 "\n"
 	 "Plugin(DisplacementRaise) is executed in-place.\n");
 }
diff --git a/Plugin/Eigenvectors.cpp b/Plugin/Eigenvectors.cpp
index 7fcc3874dd0e85881bf58b1e54724fd95a363bb3..1c2f397a353ffd77854ee71642f2226ec2ab2708 100644
--- a/Plugin/Eigenvectors.cpp
+++ b/Plugin/Eigenvectors.cpp
@@ -1,4 +1,4 @@
-// $Id: Eigenvectors.cpp,v 1.4 2005-01-01 19:35:39 geuzaine Exp $
+// $Id: Eigenvectors.cpp,v 1.5 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -58,15 +58,14 @@ void GMSH_EigenvectorsPlugin::getInfos(char *author, char *copyright, char *help
   strcpy(author, "C. Geuzaine (geuz@geuz.org)");
   strcpy(copyright, "DGR (www.multiphysics.com)");
   strcpy(help_text,
-         "Plugin(Eigenvectors) computes the three\n"
-	 "(right) eigenvectors of each tensor in the\n"
-	 "view `iView' and sorts them according to\n"
-	 "the value of the associated eigenvalues.\n"
-	 "If `ScaleByEigenvalues' is set, each\n"
-	 "eigenvector is scaled by its associated\n"
-	 "eigenvalue. The plugin gives an error if\n"
-	 "the eigenvectors are complex. If `iView'\n"
-	 "< 0, the plugin is run on the current view.\n"
+         "Plugin(Eigenvectors) computes the three (right)\n"
+	 "eigenvectors of each tensor in the view `iView'\n"
+	 "and sorts them according to the value of the\n"
+	 "associated eigenvalues. If `ScaleByEigenvalues'\n"
+	 "is set, each eigenvector is scaled by its\n"
+	 "associated eigenvalue. The plugin gives an error\n"
+	 "if the eigenvectors are complex. If `iView' < 0,\n"
+	 "the plugin is run on the current view.\n"
 	 "\n"
 	 "Plugin(Eigenvectors) creates three new\n"
 	 "vector views.\n");
diff --git a/Plugin/Levelset.cpp b/Plugin/Levelset.cpp
index ba940eadf3cff1f60ca3391950ce9a6332416af6..510df477a61890d121441d97858771e95bd248c1 100644
--- a/Plugin/Levelset.cpp
+++ b/Plugin/Levelset.cpp
@@ -1,4 +1,4 @@
-// $Id: Levelset.cpp,v 1.23 2005-01-01 19:35:39 geuzaine Exp $
+// $Id: Levelset.cpp,v 1.24 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -39,39 +39,139 @@ GMSH_LevelsetPlugin::GMSH_LevelsetPlugin()
   _valueView = -1; // use same view for levelset and field data
   _valueTimeStep = -1; // use same time step in levelset and field data views
   _recurLevel = 4;
-  _targetError = 0;
+  _targetError = 0.;
+  _extractVolume = 0; // to create isovolumes (keep all elements < or > levelset)
   _orientation = GMSH_LevelsetPlugin::NONE;
 }
 
-static void affect(double *xpi, double *ypi, double *zpi, double valpi[12][9], int i,
-		   double *xp, double *yp, double *zp, double valp[12][9], int j,
-		   int nb)
+static void affect(double *xpi, double *ypi, double *zpi, double valpi[12][9], int epi[12],
+		   int i,
+		   double *xp, double *yp, double *zp, double valp[12][9], int ep[12],
+		   int j, int nb)
 {
   xpi[i] = xp[j];
   ypi[i] = yp[j];
   zpi[i] = zp[j];
   for(int k = 0; k < nb; k++)
     valpi[i][k] = valp[j][k];
+  epi[i] = ep[j];
 }
 
-int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, 
-				      int nbNod, int nbEdg, int exn[12][2],
-				      double *x, double *y, double *z, 
-				      double *iVal, int iNbComp,
-				      double *dVal, int dNbComp,
-				      vector<Post_View *> out)
+static void removeIdenticalNodes(int *np, int nbComp, 
+				 double xp[12], double yp[12], double zp[12], 
+				 double valp[12][9], int ep[12])
 {
-  double levels[8], scalarVal[8];
+  double xpi[12], ypi[12], zpi[12], valpi[12][9];
+  int epi[12];
+
+  affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, nbComp);
+  int npi = 1;
+  for(int j = 1; j < *np; j++) {
+    for(int i = 0; i < npi; i++) {
+      if(fabs(xp[j] - xpi[i]) < 1.e-12 &&
+	 fabs(yp[j] - ypi[i]) < 1.e-12 &&
+	 fabs(zp[j] - zpi[i]) < 1.e-12) {
+	break;
+      }
+      if(i == npi-1) {
+	affect(xpi, ypi, zpi, valpi, epi, npi, xp, yp, zp, valp, ep, j, nbComp);
+	npi++;
+	break;
+      }
+    }
+  }
+  for(int i = 0; i < npi; i++)
+    affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, nbComp);
+  *np = npi;
+}
 
-  // compute the value of the levelset function at each node
+static void reorderQuad(int nbComp, 
+			double xp[12], double yp[12], double zp[12], 
+			double valp[12][9], int ep[12])
+{
+  double xpi[1], ypi[1], zpi[1], valpi[1][9];
+  int epi[12];
+  affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, nbComp);
+  affect(xp, yp, zp, valp, ep, 3, xp, yp, zp, valp, ep, 2, nbComp);
+  affect(xp, yp, zp, valp, ep, 2, xpi, ypi, zpi, valpi, epi, 0, nbComp);
+}
+
+static void reorderPrism(int nbComp, 
+			 double xp[12], double yp[12], double zp[12], 
+			 double valp[12][9], int ep[12],
+			 int nbCut, int exn[12][2])
+{
+  double xpi[6], ypi[6], zpi[6], valpi[6][9];
+  int epi[12];
+
+  if(nbCut == 3){
+    // 3 first nodes come from zero levelset intersection, next 3 are
+    // endpoints of relative edges
+    affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, nbComp);
+    affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, nbComp);
+    affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 5, nbComp);
+    for(int i = 0; i < 3; i++){
+      int edgecut = ep[i]-1;
+      for(int j = 0; j < 3; j++){
+	int p = -epi[j]-1;
+	if(exn[edgecut][0] == p || exn[edgecut][1] == p)
+	  affect(xp, yp, zp, valp, ep, 3+i, xpi, ypi, zpi, valpi, epi, j, nbComp);	  
+      }
+    }
+  }
+  else if(nbCut == 4){
+    // 4 first nodes come from zero levelset intersection
+    affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, nbComp);
+    int edgecut = ep[0]-1;
+    int p0 = -ep[4]-1;
+
+    if(exn[edgecut][0] == p0 || exn[edgecut][1] == p0){
+      affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, nbComp);
+      if(exn[ep[1]-1][0] == p0 || exn[ep[1]-1][1] == p0){
+	affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+      }
+      else{
+	affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+      }
+    }
+    else{
+      affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 5, nbComp);
+      if(exn[ep[1]-1][0] == p0 || exn[ep[1]-1][1] == p0){
+	affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+      }
+      else{
+	affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, nbComp);
+	affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+      }
+    }
+    for(int i = 0; i < 6; i++)
+      affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, nbComp);
+  }
+}
+ 
+void GMSH_LevelsetPlugin::evalLevelset(int nbNod, int nbComp,
+				       double *x, double *y, double *z, double *val,
+				       double *levels, double *scalarVal)
+{
   if(_valueIndependent) {
     for(int k = 0; k < nbNod; k++)
       levels[k] = levelset(x[k], y[k], z[k], 0.0);
   }
   else{
     for(int k = 0; k < nbNod; k++) {
-      double *vals = &iVal[iNbComp * k];
-      switch(iNbComp) {
+      double *vals = &val[nbComp * k];
+      switch(nbComp) {
       case 1: // scalar
 	scalarVal[k] = vals[0];
 	break;
@@ -85,10 +185,125 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep,
       levels[k] = levelset(x[k], y[k], z[k], scalarVal[k]);
     }
   }
+}
+
+void GMSH_LevelsetPlugin::addElement(int timeStep, int np, int nbEdg, int dNbComp,
+				     double xp[12], double yp[12], double zp[12],
+				     double valp[12][9], vector<Post_View *> out)
+{
+  // select the output view
+  Post_View *view = _valueIndependent ? out[0] : out[timeStep];
+  List_T *list;
+  int *nbPtr;
+  switch(np){
+  case 1:
+    if(dNbComp == 1)      { list = view->SP; nbPtr = &view->NbSP; }
+    else if(dNbComp == 3) { list = view->VP; nbPtr = &view->NbVP; }
+    else                  { list = view->TP; nbPtr = &view->NbTP; }
+    break;
+  case 2:
+    if(dNbComp == 1)      { list = view->SL; nbPtr = &view->NbSL; }
+    else if(dNbComp == 3) { list = view->VL; nbPtr = &view->NbVL; }
+    else                  { list = view->TL; nbPtr = &view->NbTL; }
+    break;
+  case 3:
+    if(dNbComp == 1)      { list = view->ST; nbPtr = &view->NbST; }
+    else if(dNbComp == 3) { list = view->VT; nbPtr = &view->NbVT; }
+    else                  { list = view->TT; nbPtr = &view->NbTT; }
+    break;
+  case 4:
+    if(!_extractVolume || nbEdg <= 4){
+      if(dNbComp == 1)      { list = view->SQ; nbPtr = &view->NbSQ; }
+      else if(dNbComp == 3) { list = view->VQ; nbPtr = &view->NbVQ; }
+      else                  { list = view->TQ; nbPtr = &view->NbTQ; }
+    }
+    else{
+      if(dNbComp == 1)      { list = view->SS; nbPtr = &view->NbSS; }
+      else if(dNbComp == 3) { list = view->VS; nbPtr = &view->NbVS; }
+      else                  { list = view->TS; nbPtr = &view->NbTS; }
+    }
+    break;
+  case 5:
+    if(dNbComp == 1)      { list = view->SY; nbPtr = &view->NbSY; }
+    else if(dNbComp == 3) { list = view->VY; nbPtr = &view->NbVY; }
+    else                  { list = view->TY; nbPtr = &view->NbTY; }
+    break;
+  case 6:
+    if(dNbComp == 1)      { list = view->SI; nbPtr = &view->NbSI; }
+    else if(dNbComp == 3) { list = view->VI; nbPtr = &view->NbVI; }
+    else                  { list = view->TI; nbPtr = &view->NbTI; }
+    break;
+  case 8: // should never happen
+    if(dNbComp == 1)      { list = view->SH; nbPtr = &view->NbSH; }
+    else if(dNbComp == 3) { list = view->VH; nbPtr = &view->NbVH; }
+    else                  { list = view->TH; nbPtr = &view->NbTH; }
+    break;
+  default:
+    return;
+  }
+
+  // copy the elements in the output view
+  if(!timeStep || !_valueIndependent) {
+    for(int k = 0; k < np; k++) 
+      List_Add(list, &xp[k]);
+    for(int k = 0; k < np; k++)
+      List_Add(list, &yp[k]);
+    for(int k = 0; k < np; k++)
+      List_Add(list, &zp[k]);
+    (*nbPtr)++;
+  }
+  for(int k = 0; k < np; k++)
+    for(int l = 0; l < dNbComp; l++)
+      List_Add(list, &valp[k][l]);
+}
+
+void GMSH_LevelsetPlugin::nonZeroLevelset(int timeStep, 
+					  int nbNod, int nbEdg, int exn[12][2],
+					  double *x, double *y, double *z, 
+					  double *iVal, int iNbComp,
+					  double *dVal, int dNbComp,
+					  vector<Post_View *> out)
+{
+  double levels[8], scalarVal[8];
+  
+  evalLevelset(nbNod, iNbComp, x, y, z, iVal, levels, scalarVal);
+  
+  int add = 1;
+  for(int k = 0; k < nbNod; k++){
+    if((_extractVolume < 0. && levels[k] > 0.) ||
+       (_extractVolume > 0. && levels[k] < 0.)){
+      add = 0;
+      break;
+    }
+  }
   
+  if(add){
+    double xp[12], yp[12], zp[12], valp[12][9];
+    for(int k = 0; k < nbNod; k++){
+      xp[k] = x[k];
+      yp[k] = y[k];
+      zp[k] = z[k];
+      for(int l = 0; l < dNbComp; l++)
+	valp[k][l] = dVal[dNbComp * k + l];
+    }
+    addElement(timeStep, nbNod, nbEdg, dNbComp, xp, yp, zp, valp, out);
+  }
+}
+
+int GMSH_LevelsetPlugin::zeroLevelset(int timeStep, 
+				      int nbNod, int nbEdg, int exn[12][2],
+				      double *x, double *y, double *z, 
+				      double *iVal, int iNbComp,
+				      double *dVal, int dNbComp,
+				      vector<Post_View *> out)
+{
+  double levels[8], scalarVal[8];
+
+  evalLevelset(nbNod, iNbComp, x, y, z, iVal, levels, scalarVal);
+
   // interpolate the zero levelset and the field to plot on it
   double xp[12], yp[12], zp[12], valp[12][9];
-  int np = 0;
+  int np = 0, ep[12];
   for(int k = 0; k < nbEdg; k++) {
     if(levels[exn[k][0]] * levels[exn[k][1]] <= 0.0) {
       if(iVal && dVal) {
@@ -99,6 +314,7 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep,
 	for(int l = 0; l < dNbComp; l++)
 	  valp[np][l] = coef * (val2[l] - val1[l]) + val1[l];
       }
+      ep[np] = k+1;
       np++;
     }
   }
@@ -106,48 +322,24 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep,
   if(!iVal || !dVal)
     return np;
 
-  double xpi[12], ypi[12], zpi[12], valpi[12][9];
-
   // Remove identical nodes (this can happen if an edge actually
-  // belongs to the zero levelset, i.e., if levels[] * levels[] ==
-  // 0). We should be doing this even for np < 4, but it would slow us
-  // down even more... (And we don't really care if some nodes in a
-  // postprocessing element are identical.)
-  if(np > 4) {
-    int npi;
-    affect(xpi, ypi, zpi, valpi, 0, xp, yp, zp, valp, 0, dNbComp);
-    npi = 1;
-    for(int j = 1; j < np; j++) {
-      for(int i = 0; i < npi; i++) {
-	if(fabs(xp[j] - xpi[i]) < 1.e-12 &&
-	   fabs(yp[j] - ypi[i]) < 1.e-12 &&
-	   fabs(zp[j] - zpi[i]) < 1.e-12) {
-	  break;
-	}
-	if(i == npi-1) {
-	  affect(xpi, ypi, zpi, valpi, i+1, xp, yp, zp, valp, j, dNbComp);
-	  npi++;
-	}
-      }
-    }
-    for(int i = 0; i < npi; i++)
-      affect(xp, yp, zp, valp, i, xpi, ypi, zpi, valpi, i, dNbComp);
-    np = npi;
-  }
+  // belongs to the zero levelset, i.e., if levels[] * levels[] == 0)
+  if(np > 1)
+    removeIdenticalNodes(&np, dNbComp, xp, yp, zp, valp, ep);
 
-  // can't deal with this--just return...
-  if(np < 1 || np > 4)
+  if(nbEdg > 4 && np < 3) // 3D input should only lead to 2D output
+    return 0;
+  else if(nbEdg > 1 && np < 2) // 2D input should only lead to 1D output
+    return 0;
+  else if(np < 1 || np > 4) // can't deal with this
     return 0;
 
-  // avoid ``butterflies''
-  if(np == 4) {
-    affect(xpi, ypi, zpi, valpi, 0, xp, yp, zp, valp, 3, dNbComp);
-    affect(xp, yp, zp, valp, 3, xp, yp, zp, valp, 2, dNbComp);
-    affect(xp, yp, zp, valp, 2, xpi, ypi, zpi, valpi, 0, dNbComp);
-  }
+  // avoid "butterflies"
+  if(np == 4)
+    reorderQuad(dNbComp, xp, yp, zp, valp, ep);
       
   // orient the triangles and the quads to get the normals right
-  if(np == 3 || np == 4) {
+  if(!_extractVolume && (np == 3 || np == 4)) {
     if(!timeStep || !_valueIndependent) {
       // test this only once for spatially-fixed views
       double v1[3] = { xp[2] - xp[0], yp[2] - yp[0], zp[2] - zp[0] };
@@ -173,51 +365,40 @@ int GMSH_LevelsetPlugin::zeroLevelset(int timeStep,
       }
     }
     if(_invert > 0.) {
+      double xpi[12], ypi[12], zpi[12], valpi[12][9];
+      int epi[12];
       for(int k = 0; k < np; k++)
-	affect(xpi, ypi, zpi, valpi, k, xp, yp, zp, valp, k, dNbComp);
+	affect(xpi, ypi, zpi, valpi, epi, k, xp, yp, zp, valp, ep, k, dNbComp);
       for(int k = 0; k < np; k++)
-	affect(xp, yp, zp, valp, k, xpi, ypi, zpi, valpi, np-k-1, dNbComp);
+	affect(xp, yp, zp, valp, ep, k, xpi, ypi, zpi, valpi, epi, np-k-1, dNbComp);
     }
   }
 
-  // select the output view
-  Post_View *view = _valueIndependent ? out[0] : out[timeStep];
-  List_T *list;
-  int *nbPtr;
-  if(np == 1) {
-    if(dNbComp == 1)      { list = view->SP; nbPtr = &view->NbSP; }
-    else if(dNbComp == 3) { list = view->VP; nbPtr = &view->NbVP; }
-    else                  { list = view->TP; nbPtr = &view->NbTP; }
-  }
-  else if(np == 2) {
-    if(dNbComp == 1)      { list = view->SL; nbPtr = &view->NbSL; }
-    else if(dNbComp == 3) { list = view->VL; nbPtr = &view->NbVL; }
-    else                  { list = view->TL; nbPtr = &view->NbTL; }
-  }
-  else if(np == 3) {
-    if(dNbComp == 1)      { list = view->ST; nbPtr = &view->NbST; }
-    else if(dNbComp == 3) { list = view->VT; nbPtr = &view->NbVT; }
-    else                  { list = view->TT; nbPtr = &view->NbTT; }
-  }
-  else{
-    if(dNbComp == 1)      { list = view->SQ; nbPtr = &view->NbSQ; }
-    else if(dNbComp == 3) { list = view->VQ; nbPtr = &view->NbVQ; }
-    else                  { list = view->TQ; nbPtr = &view->NbTQ; }
-  }
-  
-  // copy the elements in the output view
-  if(!timeStep || !_valueIndependent) {
-    for(int k = 0; k < np; k++) 
-      List_Add(list, &xp[k]);
-    for(int k = 0; k < np; k++)
-      List_Add(list, &yp[k]);
-    for(int k = 0; k < np; k++)
-      List_Add(list, &zp[k]);
-    (*nbPtr)++;
+  // if we compute isovolumes, add the nodes on the chosen side
+  if(_extractVolume){
+    int nbCut = np;
+    for(int k = 0; k < nbNod; k++){
+      if((_extractVolume < 0. && levels[k] < 0.0) ||
+	 (_extractVolume > 0. && levels[k] > 0.0)){
+	xp[np] = x[k];
+	yp[np] = y[k];
+	zp[np] = z[k];
+	for(int l = 0; l < dNbComp; l++)
+	  valp[np][l] = dVal[dNbComp * k + l];
+	ep[np] = -(k+1); // node num!
+	np++;
+      }
+    }
+    removeIdenticalNodes(&np, dNbComp, xp, yp, zp, valp, ep);
+    if(np == 4 && nbEdg <= 4)
+      reorderQuad(dNbComp, xp, yp, zp, valp, ep);
+    if(np == 6)
+      reorderPrism(dNbComp, xp, yp, zp, valp, ep, nbCut, exn);
+    if(np > 8) // can't deal with this
+      return 0;
   }
-  for(int k = 0; k < np; k++)
-    for(int l = 0; l < dNbComp; l++)
-      List_Add(list, &valp[k][l]);
+
+  addElement(timeStep, np, nbEdg, dNbComp, xp, yp, zp, valp, out);
   
   return 0;
 }
@@ -263,6 +444,9 @@ void GMSH_LevelsetPlugin::executeList(Post_View * iView, List_T * iList,
 						     dNbComp * nbNod * dTS);
 	  zeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z, 
 		       iVal, iNbComp, dVal, dNbComp, out);
+	  if(_extractVolume)
+	    nonZeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z, 
+			    iVal, iNbComp, dVal, dNbComp, out);
 	}
       }
     }
@@ -299,10 +483,27 @@ void GMSH_LevelsetPlugin::executeList(Post_View * iView, List_T * iList,
 		dDec.decompose(k, x, y, z, dVal, xNew, yNew, zNew, dValNew);
 		zeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri, 
 			     xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
+		if(_extractVolume)
+		  nonZeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri, 
+				  xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
 	      }
 	    }
 	  }
 	}
+	else if(_extractVolume){
+	  for(int iTS = 0; iTS < iView->NbTimeStep; iTS++) {
+	    int dTS = (dTimeStep < 0) ? iTS : dTimeStep;
+	    // don't compute the zero levelset of the value view
+	    if(dTimeStep < 0 || iView != dView || dTS != iTS) {
+	      double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbNod + 
+							 iNbComp * nbNod * iTS); 
+	      double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod + 
+							 dNbComp * nbNod * dTS);
+	      nonZeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z, 
+			      iVal, iNbComp, dVal, dNbComp, out);
+	    }
+	  }
+	}
       }
       else{
 	// since we generate one view for each time step, we can
@@ -313,21 +514,27 @@ void GMSH_LevelsetPlugin::executeList(Post_View * iView, List_T * iList,
 	  if(dTimeStep < 0 || iView != dView || dTS != iTS) {
 	    double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbNod +
 						       iNbComp * nbNod * iTS); 
+	    double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod +
+						       dNbComp * nbNod * dTS);
 	    if(zeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z, iVal, iNbComp, 
 			    NULL, 0, out)) {
-	      double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod +
-							 dNbComp * nbNod * dTS);
 	      for(int k = 0; k < iDec.numSimplices(); k++) {
 		iDec.decompose(k, x, y, z, iVal, xNew, yNew, zNew, iValNew);
 		dDec.decompose(k, x, y, z, dVal, xNew, yNew, zNew, dValNew);
 		zeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri, 
 			     xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
+		if(_extractVolume)
+		  nonZeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri, 
+				  xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
 	      }
 	    }
+	    else if(_extractVolume)
+	      nonZeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z, 
+			      iVal, iNbComp, dVal, dNbComp, out);
 	  }
 	}
       }
-
+      
       delete [] iValNew;
       delete [] dValNew;
     }
@@ -351,7 +558,6 @@ Post_View *GMSH_LevelsetPlugin::execute(Post_View * v)
   if (v->adaptive && v->NbSH)
       v->setAdaptiveResolutionLevel ( _recurLevel , this );
   
-
   if(_valueView < 0) {
     w = v;
   }
diff --git a/Plugin/Levelset.h b/Plugin/Levelset.h
index 4aa82b6cee25d36212abb05927e552a29cae30e5..c4283a4f0f4abde7d6ad3531bfc808ae58b9d34c 100644
--- a/Plugin/Levelset.h
+++ b/Plugin/Levelset.h
@@ -32,10 +32,20 @@ public:
   virtual double levelset(double x, double y, double z, double val) const = 0;
 protected:
   double _ref[3], _targetError;
-  int _valueTimeStep, _valueView, _valueIndependent, _recurLevel;  
+  int _valueTimeStep, _valueView, _valueIndependent, _recurLevel, _extractVolume;  
   ORIENTATION _orientation;
 private:
   double _invert;
+  void addElement(int timeStep, int np, int nbEdg, int dNbComp,
+		  double xp[12], double yp[12], double zp[12],
+		  double valp[12][9], vector<Post_View *> out);
+  void evalLevelset(int nbNod, int nbComp,
+		    double *x, double *y, double *z, double *val,
+		    double *levels, double *scalarVal);
+  void nonZeroLevelset(int timeStep, int nbVert, int nbEdg, int exn[12][2],
+		       double *x, double *y, double *z, 
+		       double *iVal, int iNbComp, double *dVal, int dNbComp,
+		       vector<Post_View*> out);
   int zeroLevelset(int timeStep, int nbVert, int nbEdg, int exn[12][2],
 		   double *x, double *y, double *z, 
 		   double *iVal, int iNbComp, double *dVal, int dNbComp,
diff --git a/Plugin/Makefile b/Plugin/Makefile
index afd3e34138dfb38b7bc603d2a3878ec158870c41..9ebcc7e5e643a20858fdda18beff2e91f85bb49f 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.75 2005-01-01 19:35:39 geuzaine Exp $
+# $Id: Makefile,v 1.76 2005-01-03 04:09:27 geuzaine Exp $
 #
 # Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 #
@@ -210,7 +210,7 @@ DecomposeInSimplex.o: DecomposeInSimplex.cpp Plugin.h ../Common/Options.h \
   ../DataStr/List.h ../Common/VertexArray.h ../Common/SmoothNormals.h \
   ../Common/GmshMatrix.h ../Common/AdaptiveViews.h DecomposeInSimplex.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../Common/Context.h \
-  ../DataStr/Malloc.h
+  ../DataStr/Malloc.h ../Numeric/Numeric.h
 Evaluate.o: Evaluate.cpp Plugin.h ../Common/Options.h ../Common/Message.h \
   ../Common/Views.h ../Common/ColorTable.h ../DataStr/List.h \
   ../Common/VertexArray.h ../Common/SmoothNormals.h \
diff --git a/Plugin/Remove.cpp b/Plugin/Remove.cpp
index de694c140e5642a606851fc9ace5139de675f4f0..f44b8c91c749bf254695614f2bc5e8e9a5bf5a1e 100644
--- a/Plugin/Remove.cpp
+++ b/Plugin/Remove.cpp
@@ -1,4 +1,4 @@
-// $Id: Remove.cpp,v 1.3 2005-01-01 19:35:39 geuzaine Exp $
+// $Id: Remove.cpp,v 1.4 2005-01-03 04:09:27 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -70,8 +70,8 @@ void GMSH_RemovePlugin::getInfos(char *author, char *copyright,
   strcpy(copyright, "DGR (www.multiphysics.com)");
   strcpy(help_text,
          "Plugin(Remove) removes the marked items\n"
-	 "from the view `iView'. If `iView' < 0, the\n"
-	 "plugin is run on the current view.\n"
+	 "from the view `iView'. If `iView' < 0, the plugin\n"
+	 "is run on the current view.\n"
 	 "\n"
 	 "Plugin(Remove) is executed in-place.\n");
 }
diff --git a/doc/VERSIONS b/doc/VERSIONS
index 476d2eb7a3700605263cd7afa4badafa495f58ec..ef845b641b374e157e16b201aa6d363c227fce14 100644
--- a/doc/VERSIONS
+++ b/doc/VERSIONS
@@ -1,4 +1,7 @@
-$Id: VERSIONS,v 1.295 2005-01-01 18:59:07 geuzaine Exp $
+$Id: VERSIONS,v 1.296 2005-01-03 04:09:27 geuzaine Exp $
+
+New since 1.58: All levelset-based plugins can now also compute
+isovolumes;
 
 New in 1.58: fixed UNIX socket interface on Windows (broken by the TCP
 solver patch in 1.57); bumped version number of default
diff --git a/doc/texinfo/opt_plugin.texi b/doc/texinfo/opt_plugin.texi
index a0003b9a8584f0454b352efe6f289b8f4fc8bdd6..fa3cb60eafceec74bee669b63f184915d74f2670 100644
--- a/doc/texinfo/opt_plugin.texi
+++ b/doc/texinfo/opt_plugin.texi
@@ -1,14 +1,13 @@
 @ftable @code
 @item Plugin(Annotate)
 Plugin(Annotate) adds a text string of size
-`FontSize' in the view `iView'. If `3D' is
-equal to 1, the plugin inserts the string
-in model coordinates at the position (`X',`Y',`Z').
-If `3D' is equal to 0, the plugin inserts the
-string in screen coordinates at the position
-(`X',`Y'), and aligns it according to `Align'.
-If `iView' < 0, the plugin is run on the current
-view.
+`FontSize' in the view `iView'. If `3D' is equal
+to 1, the plugin inserts the string in model
+coordinates at the position (`X',`Y',`Z'). If
+`3D' is equal to 0, the plugin inserts the string
+in screen coordinates at the position (`X',`Y'),
+and aligns it according to `Align'. If `iView'
+< 0, the plugin is run on the current view.
 
 Plugin(Annotate) is executed in-place.
 
@@ -76,14 +75,18 @@ Default value: @code{-1}
 @end table
 
 @item Plugin(CutMap)
-Plugin(CutMap) extracts the isovalue surface of
-value `A' from the view `iView' and draws the
+Plugin(CutMap) extracts the isosurface of value
+`A' from the view `iView' and draws the
 `dTimeStep'-th value of the view `dView' on this
 isovalue surface. If `iView' < 0, the plugin is
 run on the current view. If `dTimeStep' < 0, the
 plugin uses, for each time step in `iView', the
 corresponding time step in `dView'. If `dView'
 < 0, the plugin uses `iView' as the value source.
+If `extractVolume' is nonzero, the plugin
+extracts the isovolume with values smaller (if
+`extractVolume' > 0) or greater (if `extractVolume'
+< 0) than the isosurface `A'.
 
 Plugin(CutMap) creates as many views as there
 are time steps in `iView'.
@@ -96,12 +99,14 @@ Default value: @code{1}
 Default value: @code{-1}
 @item dView
 Default value: @code{-1}
-@item iView
-Default value: @code{-1}
+@item extractVolume
+Default value: @code{0}
 @item recurLevel
 Default value: @code{4}
 @item targetError
 Default value: @code{0}
+@item iView
+Default value: @code{-1}
 @end table
 
 @item Plugin(CutParametric)
@@ -141,8 +146,10 @@ Default value: @code{-1}
 @item Plugin(CutPlane)
 Plugin(CutPlane) cuts the view `iView' with
 the plane `A'*X + `B'*Y + `C'*Z + `D' = 0. If
-`iView' < 0, the plugin is run on the current
-view.
+`extractVolume' is nonzero, the plugin extracts
+the elements on one side of the plane (depending
+on the sign of `extractVolume'). If `iView' < 0,
+the plugin is run on the current view.
 
 Plugin(CutPlane) creates one new view.
 
@@ -156,19 +163,23 @@ Default value: @code{0}
 Default value: @code{0}
 @item D
 Default value: @code{-0.01}
-@item iView
-Default value: @code{-1}
+@item extractVolume
+Default value: @code{0}
 @item recurLevel
 Default value: @code{4}
 @item targetError
 Default value: @code{0}
+@item iView
+Default value: @code{-1}
 @end table
 
 @item Plugin(CutSphere)
 Plugin(CutSphere) cuts the view `iView' with the
 sphere (X-`Xc')^2 + (Y-`Yc')^2 + (Z-`Zc')^2 = `R'^2.
-If `iView' < 0, the plugin is run on the current
-view.
+If `extractVolume' is nonzero, the plugin extracts
+the elements inside (if `extractVolume' < 0) or
+outside (if `extractVolume' > 0) the sphere. If
+`iView' < 0, the plugin is run on the current view.
 
 Plugin(CutSphere) creates one new view.
 
@@ -182,10 +193,12 @@ Default value: @code{0}
 Default value: @code{0}
 @item R
 Default value: @code{0.25}
-@item iView
-Default value: @code{-1}
+@item extractVolume
+Default value: @code{0}
 @item recurLevel
 Default value: @code{4}
+@item iView
+Default value: @code{-1}
 @end table
 
 @item Plugin(DecomposeInSimplex)
@@ -209,11 +222,10 @@ Plugin(DisplacementRaise) transforms the
 coordinates of the elements in the view `iView'
 using the vectorial values (the displacements)
 stored in the `dTimeStep'-th time step of the
-view `dView'. If `iView' < 0, the plugin is
-run on the current view. If `dView' < 0, the
-plugin looks for the displacements in the
-view located just after `iView' in the view
-list.
+view `dView'. If `iView' < 0, the plugin is run
+on the current view. If `dView' < 0, the plugin
+looks for the displacements in the view located
+just after `iView' in the view list.
 
 Plugin(DisplacementRaise) is executed in-place.
 
@@ -230,15 +242,14 @@ Default value: @code{-1}
 @end table
 
 @item Plugin(Eigenvectors)
-Plugin(Eigenvectors) computes the three
-(right) eigenvectors of each tensor in the
-view `iView' and sorts them according to
-the value of the associated eigenvalues.
-If `ScaleByEigenvalues' is set, each
-eigenvector is scaled by its associated
-eigenvalue. The plugin gives an error if
-the eigenvectors are complex. If `iView'
-< 0, the plugin is run on the current view.
+Plugin(Eigenvectors) computes the three (right)
+eigenvectors of each tensor in the view `iView'
+and sorts them according to the value of the
+associated eigenvalues. If `ScaleByEigenvalues'
+is set, each eigenvector is scaled by its
+associated eigenvalue. The plugin gives an error
+if the eigenvectors are complex. If `iView' < 0,
+the plugin is run on the current view.
 
 Plugin(Eigenvectors) creates three new
 vector views.
@@ -407,8 +418,8 @@ Default value: @code{-1}
 
 @item Plugin(Remove)
 Plugin(Remove) removes the marked items
-from the view `iView'. If `iView' < 0, the
-plugin is run on the current view.
+from the view `iView'. If `iView' < 0, the plugin
+is run on the current view.
 
 Plugin(Remove) is executed in-place.