diff --git a/Box/Main.cpp b/Box/Main.cpp
index 4cac5ffefb8cd46f35b8a45fbc1f7e453f1849f1..f03db52f504f4bda0f9e145a84ad4be1bdf0fe6b 100644
--- a/Box/Main.cpp
+++ b/Box/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.34 2003-12-02 01:31:22 geuzaine Exp $
+// $Id: Main.cpp,v 1.35 2003-12-07 00:23:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -106,6 +106,7 @@ int main(int argc, char *argv[])
   M.Surfaces = NULL;
   M.Volumes = NULL;
   M.PhysicalGroups = NULL;
+  M.Partitions = NULL;
   M.Metric = NULL;
 
   signal(SIGINT, Signal);
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index c0db1eda146b5b680f9944af5480cab9fba2a488..891a6d366a43d41c63eb030ea1e11b7370d0db2d 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -696,7 +696,7 @@ StringXNumber MeshOptions_Number[] = {
   { F|O, "CharacteristicLengthFactor" , opt_mesh_lc_factor , 1.0 ,
     "Factor applied to all characteristic lengths (and background meshes)" },
   { F|O, "ColorCarousel" , opt_mesh_color_carousel , 1. ,
-    "Use a `color by region number' coloring scheme" },
+    "Use a `color by region' coloring scheme (0=no, 1=by geometrical entity, 2=by partition)" },
   { F|O, "ColorScheme" , opt_mesh_color_scheme , 0. , 
     "Default mesh color scheme (0, 1 or 2)" },
   { F|O, "ConstrainedBackgroundMesh" , opt_mesh_constrained_bgmesh, 0. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index ffdedf92748838f13493beca8df54a5fd7d910cb..3095d2e9bbd04fda5dbfa608405ce720fb0972e0 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.126 2003-12-04 02:10:30 geuzaine Exp $
+// $Id: Options.cpp,v 1.127 2003-12-07 00:23:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -3412,8 +3412,11 @@ double opt_mesh_color_carousel(OPT_ARGS_NUM)
     CTX.mesh.changed = 1;
   }
 #if defined(HAVE_FLTK)
-  if(WID && (action & GMSH_GUI))
-    WID->mesh_butt[17]->value(CTX.mesh.color_carousel);
+  if(WID && (action & GMSH_GUI)){
+    WID->mesh_butt[17]->value(CTX.mesh.color_carousel==0);
+    WID->mesh_butt[18]->value(CTX.mesh.color_carousel==1);
+    WID->mesh_butt[19]->value(CTX.mesh.color_carousel==2);
+  }
 #endif
   return CTX.mesh.color_carousel;
 }
diff --git a/Common/Visibility.cpp b/Common/Visibility.cpp
index 9758ed52dfaed697a4c0a5c9bfdbaf0d2f6bb918..7787b77381a7e86d216d6361a6e5fb4501af8f64 100644
--- a/Common/Visibility.cpp
+++ b/Common/Visibility.cpp
@@ -1,4 +1,4 @@
-// $Id: Visibility.cpp,v 1.2 2003-12-01 23:57:17 geuzaine Exp $
+// $Id: Visibility.cpp,v 1.3 2003-12-07 00:23:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -29,7 +29,7 @@
 
 extern Mesh *THEM;
 
-static List_T *ElementaryEntities = NULL, *PhysicalEntities = NULL;
+static List_T *ElementaryEntities = NULL, *PhysicalEntities = NULL, *Partitions = NULL;
 static Tree_T *VisibleThroughPhysical[4] = { NULL, NULL, NULL, NULL };
 static List_T *NumxSymb = NULL;
 static int Sort = 1;
@@ -50,6 +50,8 @@ int Entity::Num()
     return data.surface ? data.surface->Num : 0;
   case 8:
     return data.volume ? data.volume->Num : 0;
+  case 9:
+    return data.partition ? data.partition->Num : 0;
   default:
     return 0;
   }
@@ -70,6 +72,8 @@ char *Entity::Type()
   case 4:
   case 8:
     return "Volume";
+  case 9:
+    return "Partition";
   default:
     return "Unknown";
   }
@@ -108,6 +112,8 @@ int Entity::Visible()
     return data.surface ? data.surface->Visible : 0;
   case 8:
     return data.volume ? data.volume->Visible : 0;
+  case 9:
+    return data.partition ? data.partition->Visible : 0;
   default:
     return 0;
   }
@@ -214,6 +220,11 @@ void Entity::Visible(int mode)
       break;
     data.volume->Visible = mode;
     break;
+  case 9:
+    if(!data.partition)
+      break;
+    data.partition->Visible = mode;
+    break;
   }
 }
 
@@ -326,6 +337,8 @@ void Entity::RecurVisible()
         break;
       Recur(data.volume, mode);
       break;
+    case 9:
+      break;
     }
   }
 }
@@ -401,6 +414,16 @@ static void AddPhysical(void *a, void *b)
   List_Add(PhysicalEntities, &e);
 }
 
+static void AddPartition(void *a, void *b)
+{
+  MeshPartition *p = *(MeshPartition **) a;
+  Entity e;
+  e.type = 9;
+  e.data.partition = p;
+  e.str = NULL;
+  List_Add(Partitions, &e);
+}
+
 static void AddVertex(void *a, void *b)
 {
   Vertex *v = *(Vertex **) a;
@@ -467,6 +490,11 @@ List_T *GetVisibilityList(int type)
   else
     List_Reset(PhysicalEntities);
 
+  if(!Partitions)
+    Partitions = List_Create(100, 100, sizeof(Entity));
+  else
+    List_Reset(Partitions);
+
   if(!NumxSymb)
     NumxSymb = List_Create(100, 100, sizeof(NxS));
   else
@@ -475,6 +503,7 @@ List_T *GetVisibilityList(int type)
   Tree_Action(Symbol_T, addInNumxSymb);
 
   List_Action(THEM->PhysicalGroups, AddPhysical);
+  List_Action(THEM->Partitions, AddPartition);
 
   Tree_Action(THEM->Points, AddVertex);
   Tree_Action(THEM->Curves, AddCurve);
@@ -483,14 +512,16 @@ List_T *GetVisibilityList(int type)
 
   List_Sort(ElementaryEntities, CompareEntity);
   List_Sort(PhysicalEntities, CompareEntity);
+  List_Sort(Partitions, CompareEntity);
 
   switch (type) {
   case ELEMENTARY:
     return ElementaryEntities;
   case PHYSICAL:
     return PhysicalEntities;
+  case PARTITION:
   default:
-    return NULL;
+    return Partitions;
   }
 }
 
@@ -498,7 +529,13 @@ void ClearVisibilityList(int type)
 {
   int i;
   Entity *e;
-  List_T *list = (type == ELEMENTARY) ? ElementaryEntities : PhysicalEntities;
+  List_T *list;
+  switch(type){
+  case ELEMENTARY : list = ElementaryEntities; break;
+  case PHYSICAL : list = PhysicalEntities; break;
+  case PARTITION : list = Partitions; break;
+  default: return;
+  }
   for(i = 0; i < List_Nbr(list); i++) {
     e = (Entity *) List_Pointer(list, i);
     e->Visible(0);
@@ -648,6 +685,10 @@ static void vis_pyr(void *a, void *b)
 {
   (*(Pyramid **) a)->Visible = vmode;
 }
+static void vis_part(void *a, void *b)
+{
+  (*(MeshPartition **) a)->Visible = vmode;
+}
 static void vis_cur(void *a, void *b)
 {
   (*(Curve **) a)->Visible = vmode;
diff --git a/Common/Visibility.h b/Common/Visibility.h
index b90d51f294072b99c163240a7199ab28993c48f6..49f1f81db43e7cba99bc874291d763003cd21382 100644
--- a/Common/Visibility.h
+++ b/Common/Visibility.h
@@ -37,6 +37,7 @@ public :
     Surface *surface;
     Volume *volume;
     PhysicalGroup *physical;
+    MeshPartition *partition;
   } data;
   char *str;
 
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 58a66f91a88b515c09e22ecad4ee13eb02c62bbf..1422f978f96ee0664c30c7791875ddf17bdafe32 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.195 2003-12-04 02:10:31 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.196 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -1038,7 +1038,9 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
   opt_mesh_aspect(0, GMSH_SET,
                   WID->mesh_butt[14]->value()? 0 :
                   WID->mesh_butt[15]->value()? 1 : 2);
-  opt_mesh_color_carousel(0, GMSH_SET, WID->mesh_butt[17]->value());
+  opt_mesh_color_carousel(0, GMSH_SET,
+			  WID->mesh_butt[17]->value()? 0 :
+			  WID->mesh_butt[18]->value()? 1 : 2);
 
   opt_mesh_nb_smoothing(0, GMSH_SET, WID->mesh_value[0]->value());
   opt_mesh_scaling_factor(0, GMSH_SET, WID->mesh_value[1]->value());
@@ -1233,9 +1235,12 @@ void visibility_cb(CALLBACK_ARGS)
   case 0:
     type = ELEMENTARY;
     break;
-  default:
+  case 1:
     type = PHYSICAL;
     break;
+  default :
+    type = PARTITION;
+    break;
   }
   switch (WID->vis_browser_mode->value()) {
   case 0:
@@ -1269,9 +1274,12 @@ void visibility_ok_cb(CALLBACK_ARGS)
   case 0:
     ClearVisibilityList(PHYSICAL);
     break;
-  default:
+  case 1:
     ClearVisibilityList(ELEMENTARY);
     break;
+  default:
+    // partitions: do nothing
+    break;
   }
   switch (WID->vis_browser_mode->value()) {
   case 0:
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index ff4fcf228f542d3158cfefb2d037a497ee79654d..b99b63ccb7e1a650328407a6ab811cfa05bbc470 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.264 2003-12-04 02:10:31 geuzaine Exp $
+// $Id: GUI.cpp,v 1.265 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -1824,22 +1824,26 @@ void GUI::create_option_window()
     {
       Fl_Group *o = new Fl_Group(WB, WB + BH, width - 2 * WB, height - 2 * WB - BH, "Colors");
       o->hide();
-      mesh_butt[17] = new Fl_Check_Button(2 * WB, 2 * WB + 1 * BH, BW, BH, "Switch color by entity");
-      mesh_butt[17]->type(FL_TOGGLE_BUTTON);
-      mesh_butt[17]->down_box(TOGGLE_BOX);
-      mesh_butt[17]->selection_color(TOGGLE_COLOR);
+      mesh_butt[17] = new Fl_Check_Button(2 * WB, 2 * WB + 1 * BH, BW, BH, "Color by element type");
+      mesh_butt[18] = new Fl_Check_Button(2 * WB, 2 * WB + 2 * BH, BW, BH, "Color by geometrical entity");
+      mesh_butt[19] = new Fl_Check_Button(2 * WB, 2 * WB + 3 * BH, BW, BH, "Color by partition");
+      for(i = 17; i < 20; i++) {
+        mesh_butt[i]->type(FL_RADIO_BUTTON);
+        mesh_butt[i]->down_box(RADIO_BOX);
+        mesh_butt[i]->selection_color(RADIO_COLOR);
+      }
 
-      mesh_value[12] = new Fl_Value_Input(2 * WB, 2 * WB + 2 * BH, IW, BH, "Predefined color scheme");
+      mesh_value[12] = new Fl_Value_Input(2 * WB, 2 * WB + 4 * BH, IW, BH, "Predefined color scheme");
       mesh_value[12]->minimum(0);
       mesh_value[12]->maximum(2);
       mesh_value[12]->step(1);
       mesh_value[12]->align(FL_ALIGN_RIGHT);
       mesh_value[12]->callback(mesh_options_color_scheme_cb);
 
-      Fl_Scroll *s = new Fl_Scroll(2 * WB, 3 * WB + 3 * BH, IW + 20, height - 5 * WB - 3 * BH);
+      Fl_Scroll *s = new Fl_Scroll(2 * WB, 3 * WB + 5 * BH, IW + 20, height - 5 * WB - 5 * BH);
       i = 0;
       while(MeshOptions_Color[i].str) {
-        mesh_col[i] = new Fl_Button(2 * WB, 3 * WB + (3 + i) * BH, IW, BH, MeshOptions_Color[i].str);
+        mesh_col[i] = new Fl_Button(2 * WB, 3 * WB + (5 + i) * BH, IW, BH, MeshOptions_Color[i].str);
         mesh_col[i]->callback(color_cb, (void *)MeshOptions_Color[i].function);
         i++;
       }
@@ -2823,6 +2827,7 @@ void GUI::create_visibility_window()
   static Fl_Menu_Item type_table[] = {
     {"Elementary", 0, (Fl_Callback *) visibility_cb},
     {"Physical", 0, (Fl_Callback *) visibility_cb},
+    {"Partitions", 0, (Fl_Callback *) visibility_cb},
     {0}
   };
   static Fl_Menu_Item browser_mode_table[] = {
diff --git a/Fltk/Main.cpp b/Fltk/Main.cpp
index 21a80eb352da68b7763ca39df216d7ea51ded3ac..16a24f546747f20e184c12d12c98faea74a8d897 100644
--- a/Fltk/Main.cpp
+++ b/Fltk/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.53 2003-11-29 19:29:26 geuzaine Exp $
+// $Id: Main.cpp,v 1.54 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -114,6 +114,7 @@ int main(int argc, char *argv[])
   M.Surfaces = NULL;
   M.Volumes = NULL;
   M.PhysicalGroups = NULL;
+  M.Partitions = NULL;
   M.Metric = NULL;
 
   // Signal handling
diff --git a/Geo/Geo.h b/Geo/Geo.h
index c332e6cedf9df61dcb7033bb4f4aef23a0ba9a2f..d04154b1fca41a9d3413522df01ba192d57da2bf 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -24,6 +24,7 @@
 
 #define ELEMENTARY 1
 #define PHYSICAL   2
+#define PARTITION  3
 
 #define INFILE     1
 #define INSTRING   2
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index c5c42b86755e41457764223fcc059a6247b683c4..6712687824530551f4389344d5869bbd6a7010cb 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.62 2003-06-14 04:37:42 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.63 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -59,12 +59,11 @@ void draw_polygon_2d(double r, double g, double b, int n,
 
 }
 
-static int iColor;
 static DrawingColor theColor;
 
 void ColorSwitch(int i)
 {
-  glColor4ubv((GLubyte *) & CTX.color.mesh.carousel[i % 10]);
+  glColor4ubv((GLubyte *) & CTX.color.mesh.carousel[abs(i % 10)]);
 }
 
 void Draw_Mesh(Mesh * M)
@@ -123,8 +122,6 @@ void Draw_Mesh(Mesh * M)
     glEnd();
   }
 
-  iColor = 0;
-
   if(CTX.mesh.draw && CTX.render_mode != GMSH_SELECT) {
 
     static int first = 1, listnum;
@@ -197,7 +194,6 @@ void Draw_Mesh_Volumes(void *a, void *b)
 {
   Volume *v;
   v = *(Volume **) a;
-  iColor++;
   theColor = v->Color;
   // FIXME: this is the correct method, but will only work when a
   // coherent datastruct exists for volumes
@@ -222,7 +218,6 @@ void Draw_Mesh_Surfaces(void *a, void *b)
 {
   Surface *s;
   s = *(Surface **) a;
-  iColor++;
   theColor = s->Color;
   if(!(s->Visible & VIS_MESH))
     return;
@@ -250,7 +245,6 @@ void Draw_Mesh_Curves(void *a, void *b)
   c = *(Curve **) a;
   if(c->Num < 0)
     return;
-  iColor++;
   theColor = c->Color;
   if(!(c->Visible & VIS_MESH))
     return;
@@ -317,6 +311,10 @@ void Draw_Simplex_Volume(void *a, void *b)
   if(!s->V[3] || !(s->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, s->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   // FIXME: remove as soon as a coherent structure exists for volumes
   Volume *V;
   if((V = FindVolume(s->iEnt, THEM)) && !(V->Visible & VIS_MESH))
@@ -349,8 +347,9 @@ void Draw_Simplex_Volume(void *a, void *b)
     fulldraw = 1;
   }
 
-  if(CTX.mesh.color_carousel && !fulldraw)
-    ColorSwitch(s->iEnt);
+  if(CTX.mesh.color_carousel && !fulldraw){
+    ColorSwitch((CTX.mesh.color_carousel == 2) ? s->iPart : s->iEnt);
+  }
   else if(fulldraw)
     glColor4ubv((GLubyte *) & CTX.color.mesh.line);
   else
@@ -430,7 +429,7 @@ void Draw_Simplex_Volume(void *a, void *b)
   double n[4], x1x0, y1y0, z1z0, x2x0, y2y0, z2z0;
 
   if(CTX.mesh.color_carousel)
-    ColorSwitch(s->iEnt);
+    ColorSwitch((CTX.mesh.color_carousel == 2) ? s->iPart : s->iEnt);
   else
     glColor4ubv((GLubyte *) & CTX.color.mesh.tetrahedron);
 
@@ -546,7 +545,7 @@ void Draw_Simplex_Surface_Common(Simplex * s, double *pX, double *pY,
       if(theColor.type)
         glColor4ubv((GLubyte *) & theColor.mesh);
       else
-        ColorSwitch(iColor);
+        ColorSwitch((CTX.mesh.color_carousel == 2) ? s->iPart : s->iEnt);
     }
     else {
       glColor4ubv((GLubyte *) & CTX.color.mesh.line);
@@ -589,7 +588,7 @@ void Draw_Simplex_Surface_Common(Simplex * s, double *pX, double *pY,
     if(theColor.type)
       glColor4ubv((GLubyte *) & theColor.mesh);
     else
-      ColorSwitch(iColor);
+      ColorSwitch((CTX.mesh.color_carousel == 2) ? s->iPart : s->iEnt);
   }
   else {
     if(K == 3)
@@ -663,6 +662,10 @@ void Draw_Simplex_Surface_Simple(void *a, void *b)
   if(!s->V[2] || !(s->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, s->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   Draw_Simplex_Surface_Common(s, NULL, NULL, NULL, n);
 }
 
@@ -679,6 +682,10 @@ void Draw_Simplex_Surface(void *a, void *b)
   if(!s->V[2] || !(s->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, s->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   L = (s->VSUP) ? 1 : 0;
   K = (s->V[3]) ? 4 : 3;
 
@@ -769,6 +776,10 @@ void Draw_Simplex_Curves(void *a, void *b)
   if(!(s->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, s->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   Xc = 0.5 * (s->V[0]->Pos.X + s->V[1]->Pos.X);
   Yc = 0.5 * (s->V[0]->Pos.Y + s->V[1]->Pos.Y);
   Zc = 0.5 * (s->V[0]->Pos.Z + s->V[1]->Pos.Z);
@@ -795,7 +806,7 @@ void Draw_Simplex_Curves(void *a, void *b)
     if(theColor.type)
       glColor4ubv((GLubyte *) & theColor.mesh);
     else
-      ColorSwitch(iColor);
+      ColorSwitch((CTX.mesh.color_carousel == 2) ? s->iPart : s->iEnt);
   }
   else
     glColor4ubv((GLubyte *) & CTX.color.mesh.line);
@@ -849,6 +860,10 @@ void Draw_Hexahedron_Volume(void *a, void *b)
   if(!(h->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, h->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   // FIXME: remove as soon as a coherent structure exists for volumes
   Volume *V;
   if((V = FindVolume(h->iEnt, THEM)) && !(V->Visible & VIS_MESH))
@@ -869,7 +884,7 @@ void Draw_Hexahedron_Volume(void *a, void *b)
   }
 
   if(CTX.mesh.color_carousel)
-    ColorSwitch(h->iEnt);
+    ColorSwitch((CTX.mesh.color_carousel == 2) ? h->iPart : h->iEnt);
   else
     glColor4ubv((GLubyte *) & CTX.color.mesh.hexahedron);
 
@@ -990,6 +1005,10 @@ void Draw_Prism_Volume(void *a, void *b)
   if(!(p->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, p->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   // FIXME: remove as soon as a coherent structure exists for volumes
   Volume *V;
   if((V = FindVolume(p->iEnt, THEM)) && !(V->Visible & VIS_MESH))
@@ -1010,7 +1029,7 @@ void Draw_Prism_Volume(void *a, void *b)
   }
 
   if(CTX.mesh.color_carousel)
-    ColorSwitch(p->iEnt);
+    ColorSwitch((CTX.mesh.color_carousel == 2) ? p->iPart : p->iEnt);
   else
     glColor4ubv((GLubyte *) & CTX.color.mesh.prism);
 
@@ -1112,6 +1131,10 @@ void Draw_Pyramid_Volume(void *a, void *b)
   if(!(p->Visible & VIS_MESH))
     return;
 
+  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, p->iPart);
+  if(part && !(*part)->Visible)
+    return;
+
   // FIXME: remove as soon as a coherent structure exists for volumes
   Volume *V;
   if((V = FindVolume(p->iEnt, THEM)) && !(V->Visible & VIS_MESH))
@@ -1132,7 +1155,7 @@ void Draw_Pyramid_Volume(void *a, void *b)
   }
 
   if(CTX.mesh.color_carousel)
-    ColorSwitch(p->iEnt);
+    ColorSwitch((CTX.mesh.color_carousel == 2) ? p->iPart : p->iEnt);
   else
     glColor4ubv((GLubyte *) & CTX.color.mesh.pyramid);
 
diff --git a/Mesh/Create.cpp b/Mesh/Create.cpp
index 137ce51eae8741fc7dee5b64c45768290507769b..1e50f8846e24573436abb2cb5e34659fdd2c1885 100644
--- a/Mesh/Create.cpp
+++ b/Mesh/Create.cpp
@@ -1,4 +1,4 @@
-// $Id: Create.cpp,v 1.45 2003-05-29 14:36:56 geuzaine Exp $
+// $Id: Create.cpp,v 1.46 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -202,6 +202,59 @@ void Add_PhysicalGroup(int Num, int typ, List_T * intlist, Mesh * M)
   List_Add(M->PhysicalGroups, &p);
 }
 
+void Free_PhysicalGroup(void *a, void *b)
+{
+  PhysicalGroup *p = *(PhysicalGroup **) a;
+  if(p) {
+    List_Delete(p->Entities);
+    Free(p);
+    p = NULL;
+  }
+}
+
+int Add_MeshPartition(int Num, Mesh * M)
+{
+  MeshPartition P, *p, **pp;
+  p = &P;
+  p->Num = Num;
+  if((pp = (MeshPartition**)List_PQuery(M->Partitions, &p, compareMeshPartitionNum))){
+    return (*pp)->Index;
+  }
+  else{
+    p = (MeshPartition*)Malloc(sizeof(MeshPartition));
+    p->Num = Num;
+    p->Visible = VIS_GEOM | VIS_MESH;
+    p->Index = List_Nbr(M->Partitions);
+    List_Add(M->Partitions, &p);
+    return p->Index;
+  }
+}
+
+void Free_MeshPartition(void *a, void *b)
+{
+  MeshPartition *p = *(MeshPartition **) a;
+  if(p) {
+    Free(p);
+    p = NULL;
+  }
+}
+
+int compareMeshPartitionNum(const void *a, const void *b)
+{
+  MeshPartition *q, *w;
+  q = *(MeshPartition **) a;
+  w = *(MeshPartition **) b;
+  return (q->Num - w->Num);
+}
+
+int compareMeshPartitionIndex(const void *a, const void *b)
+{
+  MeshPartition *q, *w;
+  q = *(MeshPartition **) a;
+  w = *(MeshPartition **) b;
+  return (q->Index - w->Index);
+}
+
 void Add_EdgeLoop(int Num, List_T * intlist, Mesh * M)
 {
   EdgeLoop *pEL;
@@ -692,6 +745,7 @@ Hexahedron *Create_Hexahedron(Vertex * v1, Vertex * v2, Vertex * v3,
 
   h = (Hexahedron *) Malloc(sizeof(Hexahedron));
   h->iEnt = -1;
+  h->iPart = -1;
   h->Num = ++THEM->MaxSimplexNum;
   h->Visible = VIS_MESH;
   h->V[0] = v1;
@@ -723,6 +777,7 @@ Prism *Create_Prism(Vertex * v1, Vertex * v2, Vertex * v3,
 
   p = (Prism *) Malloc(sizeof(Prism));
   p->iEnt = -1;
+  p->iPart = -1;
   p->Num = ++THEM->MaxSimplexNum;
   p->Visible = VIS_MESH;
   p->V[0] = v1;
@@ -752,6 +807,7 @@ Pyramid *Create_Pyramid(Vertex * v1, Vertex * v2, Vertex * v3,
 
   p = (Pyramid *) Malloc(sizeof(Pyramid));
   p->iEnt = -1;
+  p->iPart = -1;
   p->Num = ++THEM->MaxSimplexNum;
   p->Visible = VIS_MESH;
   p->V[0] = v1;
diff --git a/Mesh/Create.h b/Mesh/Create.h
index b5db504342695efd031b73854a876e431f4296a5..f61f63d388351c63f54ff9f7d12d5ab726b4b7b8 100644
--- a/Mesh/Create.h
+++ b/Mesh/Create.h
@@ -37,21 +37,28 @@ int compareAttractor (const void *a, const void *b);
 int compareSurface (const void *a, const void *b);
 int compareVolume (const void *a, const void *b);
 int compareSxF (const void *a, const void *b);
+int compareMeshPartitionNum(const void *a, const void *b);
+int compareMeshPartitionIndex(const void *a, const void *b);
 
 Attractor * Create_Attractor (int Num, double lc1, double lc2, double Radius,
                               Vertex * v, Curve * c, Surface * s);
+
+void Add_EdgeLoop (int Num, List_T * intlist, Mesh * M);
 void Add_SurfaceLoop (int Num, List_T * intlist, Mesh * M);
+
 void Add_PhysicalGroup (int Num, int typ, List_T * intlist, Mesh * M);
-void Add_EdgeLoop (int Num, List_T * intlist, Mesh * M);
+void Free_PhysicalGroup(void *a, void *b);
 
-void End_Curve (Curve * c);
-void End_Surface (Surface * s);
+int  Add_MeshPartition(int Num, Mesh * M);
+void Free_MeshPartition(void *a, void *b);
 
 Curve *Create_Curve (int Num, int Typ, int Order, List_T * Liste,
                      List_T * Knots, int p1, int p2, double u1, double u2);
+void End_Curve (Curve * c);
 void Free_Curve(void *a, void *b);
 
 Surface * Create_Surface (int Num, int Typ);
+void End_Surface (Surface * s);
 void Free_Surface(void *a, void *b);
 
 Volume * Create_Volume (int Num, int Typ);
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 642e32ef54fb5e13cb67d81767b601aad8d712ad..47052c8a166e81549a8e2f4cf4ec6dbfbd2c694c 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.45 2003-06-14 04:37:42 geuzaine Exp $
+// $Id: Generator.cpp,v 1.46 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -212,9 +212,12 @@ void Init_Mesh(Mesh * M, int all)
   Tree_Action(M->Volumes, Free_Volume);
   Tree_Delete(M->Volumes);
 
-  //Tree_Action (M->PhysicalGroups, Free_PhysicalGroup); // todo
+  List_Action(M->PhysicalGroups, Free_PhysicalGroup);
   List_Delete(M->PhysicalGroups);
 
+  List_Action(M->Partitions, Free_MeshPartition);
+  List_Delete(M->Partitions);
+
   if(M->Metric) {
     delete M->Metric;
   }
@@ -228,6 +231,7 @@ void Init_Mesh(Mesh * M, int all)
   M->Surfaces = Tree_Create(sizeof(Surface *), compareSurface);
   M->Volumes = Tree_Create(sizeof(Volume *), compareVolume);
   M->PhysicalGroups = List_Create(5, 5, sizeof(PhysicalGroup *));
+  M->Partitions = List_Create(5, 5, sizeof(MeshPartition *));
   M->Metric = new GMSHMetric;
   M->BGM.bgm = NULL;
 
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index 8a729fa9c138f9f12b4b03c5827207433d2765b0..b7e7276b0b05068be090713ccdbc2928f2d4bb7e 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -99,11 +99,11 @@ typedef struct _DELAUNAY Delaunay, *delpeek;
 typedef int PointNumero;
 
 struct _DOC{
-  PointRecord *points;  /* points a trianguler */
-  List_T *hotpoints;    /* hotpoints */
-  int numPoints;        /* nombre de points */
-  int numTriangles;     /* nombre de triangles */
-  Delaunay *delaunay;   /* resultats 2D */
+  PointRecord *points;  // points to triangulate
+  List_T *hotpoints;    // hotpoints
+  int numPoints;        // number of points
+  int numTriangles;     // number of triangles
+  Delaunay *delaunay;   // 2D results
 };
 
 typedef struct{
@@ -124,7 +124,7 @@ typedef struct{
 }IntPoint;
 
 struct _CDLIST{
-  PointNumero point_num; /* numero du point */
+  PointNumero point_num;
   DListPeek next, prev;
 };
 
@@ -197,27 +197,30 @@ class NXE{
 };
 
 typedef struct{
-  int Num;              /* Numero                                       */
-  int iEnt;             /* Entite geometrique                           */
-  char Visible;         /* Visualization flag                           */
-  Vertex *V[8];         /* 8 noeuds                                     */
-  Vertex **VSUP;        /* noeuds supplem pour les elts de degre eleves */
+  int Num; 
+  int iEnt; // parent geometrical entity
+  int iPart; // mesh partition index
+  char Visible;
+  Vertex *V[8];
+  Vertex **VSUP;
 }Hexahedron;
 
 typedef struct{
-  int Num;              /* Numero                                       */
-  int iEnt;             /* Entite geometrique                           */
-  char Visible;         /* Visualization flag                           */
-  Vertex *V[6];         /* 6 noeuds                                     */
-  Vertex **VSUP;        /* noeuds supplem pour les elts de degre eleves */
+  int Num;
+  int iEnt; // parent geometrical entity
+  int iPart; // mesh partition index
+  char Visible;
+  Vertex *V[6];
+  Vertex **VSUP;
 }Prism;
 
 typedef struct{
-  int Num;              /* Numero                                       */
-  int iEnt;             /* Entite geometrique                           */
-  char Visible;         /* Visualization flag                           */
-  Vertex *V[5];         /* 5 noeuds                                     */
-  Vertex **VSUP;        /* noeuds supplem pour les elts de degre eleves */
+  int Num;
+  int iEnt; // parent geometrical entity
+  int iPart; // mesh partition index
+  char Visible;
+  Vertex *V[5];
+  Vertex **VSUP;
 }Pyramid;
 
 typedef struct{
@@ -254,8 +257,8 @@ struct _Surf{
   double RecombineAngle;
   int ipar[5];
   int Nu, Nv;
-  List_T *Generatrices; /* Surface reglee    */
-  List_T *Control_Points; /* Patchs bicubiques */
+  List_T *Generatrices;
+  List_T *Control_Points;
   double plan[3][3];
   double invplan[3][3];
   double a, b, c, d;
@@ -270,14 +273,13 @@ struct _Surf{
   float *ku, *kv, *cp;
   struct _Surf *Support;
   CylParam Cyl;
-  Grid_T Grid;          /* Grille de recherches rapides */
+  Grid_T Grid;  // fast search grid
   ExtrudeParams *Extrude;
   STL_Data *STL; // stl representation of the surface
   POLY_rep *thePolyRep;
-  int Dirty; //flag to prevent any meshing
+  int Dirty; // flag to prevent any meshing
   DrawingColor Color;
-  /// a pointer to a solid model entity
-  void * aSolidModelEntity;
+  void * aSolidModelEntity; // pointer to a solid model entity
 };
 
 typedef struct _Surf Surface;
@@ -299,6 +301,12 @@ typedef struct{
   List_T *Entities;
 }PhysicalGroup;
 
+typedef struct{
+  int Index;
+  int Num;
+  char Visible;
+}MeshPartition;
+
 typedef struct{
   Face F;
   Face Sorted;
@@ -318,7 +326,7 @@ typedef struct {
   Tree_T *Edges;
   Tree_T *Faces;
   Tree_T *Simplexes;
-  Tree_T *Simp_Surf;//for old extrusion mesh generator
+  Tree_T *Simp_Surf; // for old extrusion mesh generator
   Tree_T *Hexahedra;
   Tree_T *Prisms;
   Tree_T *Pyramids;
@@ -331,23 +339,23 @@ typedef struct {
   int iFac;
 }exf_T;
 
-/* Structure intersection arete - Simplexe */
+// Edge-Simplex intersections
 
 typedef struct{
-  int NbIntersect;      /* nombre total d'intersections                   */
-  Edge *e;              /* arete                                          */
-  Simplex *s;           /* simplexe                                       */
-  Face *f;              /* face                                           */
-  int NbVertex;         /* nombre de noeuds du simplexe que coupe l'arete */
-  Vertex *V[12];        /* noeuds du simplexe que coupe l'arete           */
-  int iV[12];           /* noeuds du simplexe que coupe l'arete           */
-  int NbEdge;           /* nombre d'intersections arete-arete             */
-  int E[12];            /* aretes                                         */
-  Vertex *VE[12];       /* noeuds d'intersection                          */
-  int NbFace;           /* nombre d'intersections face-arete              */
-  Face *F[12];          /* faces                                          */
-  int iF[12];           /* faces                                          */
-  Vertex *VF[12];       /* position des points d'intersections face-arete */
+  int NbIntersect;      // nombre total d'intersections
+  Edge *e;              // arete
+  Simplex *s;           // simplexe
+  Face *f;              // face
+  int NbVertex;         // nombre de noeuds du simplexe que coupe l'arete
+  Vertex *V[12];        // noeuds du simplexe que coupe l'arete
+  int iV[12];           // noeuds du simplexe que coupe l'arete
+  int NbEdge;           // nombre d'intersections arete-arete
+  int E[12];            // aretes
+  Vertex *VE[12];       // noeuds d'intersection
+  int NbFace;           // nombre d'intersections face-arete
+  Face *F[12];          // faces
+  int iF[12];           // faces
+  Vertex *VF[12];       // position des points d'intersections face-arete
 }Intersection;
 
 typedef struct _Mesh Mesh;
@@ -386,7 +394,7 @@ typedef struct{
   int degre;
   CircParam Circle;
   char functu[256], functv[256], functw[256];
-  int Dirty; //flag to prevent any meshing
+  int Dirty; // flag to prevent any meshing
   DrawingColor Color;
 }Curve;
 
@@ -421,22 +429,23 @@ class MeshParameters{
 };
 
 struct _Mesh{
-  char name[256];               /* Nom du probleme                       */
-  int status;                   /* Etat actuel du maillage               */
-  Tree_T *Points;               /* Points de controle                    */
-  Tree_T *Vertices;             /* Noeuds du maillage                    */
-  Tree_T *Simplexes;            /* Simplexes                             */
-  Tree_T *Curves;               /* Courbes                               */
-  Tree_T *Surfaces;             /* Surfaces                              */
-  Tree_T *Volumes;              /* Volumes                               */
-  Tree_T *SurfaceLoops;         /* Surface Loops                         */
-  Tree_T *EdgeLoops;            /* Edge Loops                            */
-  List_T *PhysicalGroups;       /* Physical Groups                       */
-  Grid_T Grid;                  /* Grille de recherches rapides          */
-  LcField BGM;                  /* Background mesh                       */
-  double Statistics[50];        /* Mesh statistics                       */
-  int Histogram[3][NB_HISTOGRAM]; /* Quality histograms                 */
-  GMSHMetric *Metric;           /* Metric                                */
+  char name[256];
+  int status; // current state of the mesh
+  Tree_T *Points;
+  Tree_T *Vertices;
+  Tree_T *Simplexes;
+  Tree_T *Curves;
+  Tree_T *Surfaces;
+  Tree_T *Volumes;
+  Tree_T *SurfaceLoops;
+  Tree_T *EdgeLoops;
+  List_T *PhysicalGroups;
+  List_T *Partitions;
+  Grid_T Grid; // fast search grid
+  LcField BGM; // background mesh
+  double Statistics[50]; // mesh statistics
+  int Histogram[3][NB_HISTOGRAM]; // quality histograms
+  GMSHMetric *Metric;
   MeshParameters MeshParams;
   int MaxPointNum, MaxLineNum, MaxLineLoopNum, MaxSurfaceNum;
   int MaxSurfaceLoopNum, MaxVolumeNum, MaxPhysicalNum;
@@ -455,7 +464,7 @@ struct Map{
 };
 
 
-/* public functions */
+// public functions
 
 void mai3d (Mesh * M, int Asked);
 
diff --git a/Mesh/Read_Mesh.cpp b/Mesh/Read_Mesh.cpp
index fd3c802205ca9cb89760dabe0b01f91063789b31..50af7084ecbcf8f3a98a6097882bb765f78d2598 100644
--- a/Mesh/Read_Mesh.cpp
+++ b/Mesh/Read_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Read_Mesh.cpp,v 1.59 2003-11-27 02:33:31 geuzaine Exp $
+// $Id: Read_Mesh.cpp,v 1.60 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -172,7 +172,6 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 
         fscanf(fp, "%d %d %d %d %d",
 	       &Num, &Type, &Physical, &Elementary, &Nbr_Nodes);
-	//&Num, &Type, &Elementary, &Physical, &Nbr_Nodes) ;
 
         for(j = 0; j < Nbr_Nodes; j++)
           fscanf(fp, "%d", &verts[j].Num);
@@ -266,6 +265,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           simp = Create_Simplex(vertsp[0], vertsp[1], NULL, NULL);
           simp->Num = Num;
           simp->iEnt = Elementary;
+          simp->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == LGN2){
 	    simp->VSUP = (Vertex **) Malloc(1 * sizeof(Vertex *));
 	    simp->VSUP[0] = vertsp[2];
@@ -281,6 +281,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           simp = Create_Simplex(vertsp[0], vertsp[1], vertsp[2], NULL);
           simp->Num = Num;
           simp->iEnt = Elementary;
+          simp->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == TRI2){
 	    simp->VSUP = (Vertex **) Malloc(3 * sizeof(Vertex *));
 	    for(i = 0; i < 3; i++){
@@ -300,6 +301,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           simp = Create_Quadrangle(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
           simp->Num = Num;
           simp->iEnt = Elementary;
+          simp->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == QUA2){
 	    simp->VSUP = (Vertex **) Malloc(4 * sizeof(Vertex *));
 	    for(i = 0; i < 4; i++){
@@ -320,6 +322,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
           simp = Create_Simplex(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
           simp->Num = Num;
           simp->iEnt = Elementary;
+          simp->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == TET2){
 	    simp->VSUP = (Vertex **) Malloc(6 * sizeof(Vertex *));
 	    for(i = 0; i < 6; i++){
@@ -340,6 +343,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
                                   vertsp[4], vertsp[5], vertsp[6], vertsp[7]);
           hex->Num = Num;
           hex->iEnt = Elementary;
+          hex->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == HEX2){
 	    hex->VSUP = (Vertex **) Malloc(12 * sizeof(Vertex *));
 	    for(i = 0; i < 12; i++){
@@ -360,6 +364,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
                              vertsp[3], vertsp[4], vertsp[5]);
           pri->Num = Num;
           pri->iEnt = Elementary;
+          pri->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == PRI2){
 	    pri->VSUP = (Vertex **) Malloc(9 * sizeof(Vertex *));
 	    for(i = 0; i < 9; i++){
@@ -380,6 +385,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
                                vertsp[3], vertsp[4]);
           pyr->Num = Num;
           pyr->iEnt = Elementary;
+          pyr->iPart = Add_MeshPartition(Physical, M);
 	  if(Type == PYR2){
 	    pyr->VSUP = (Vertex **) Malloc(8 * sizeof(Vertex *));
 	    for(i = 0; i < 8; i++){
@@ -427,20 +433,26 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 
   if(Tree_Nbr(M->Volumes)) {
     M->status = 3;
-    M->Statistics[6] = Tree_Nbr(M->Vertices);   //incorrect, mais...
+    M->Statistics[6] = Tree_Nbr(M->Vertices);   // wrong, but...
   }
   else if(Tree_Nbr(M->Surfaces)) {
     M->status = 2;
-    M->Statistics[5] = Tree_Nbr(M->Vertices);   //incorrect, mais...
+    M->Statistics[5] = Tree_Nbr(M->Vertices);   // wrong, but...
   }
   else if(Tree_Nbr(M->Curves)) {
     M->status = 1;
-    M->Statistics[4] = Tree_Nbr(M->Vertices);   //incorrect, mais...
+    M->Statistics[4] = Tree_Nbr(M->Vertices);   // wrong, but...
   }
   else if(Tree_Nbr(M->Points))
     M->status = 0;
   else
     M->status = -1;
+
+  // For efficiency reasons, we store the partition index (and not the
+  // partition number) in the various mesh elements. We need to
+  // re-sort the list according to these indices to allow direct
+  // access through Llist_Pointer & co.
+  List_Sort(M->Partitions, compareMeshPartitionIndex);
 }
 
 
diff --git a/Mesh/Simplex.cpp b/Mesh/Simplex.cpp
index 668ded107f970947151d51545d6066f790af5ee2..37d72412ac4c8ce8dc1652cb77611877b4d9a254 100644
--- a/Mesh/Simplex.cpp
+++ b/Mesh/Simplex.cpp
@@ -1,4 +1,4 @@
-// $Id: Simplex.cpp,v 1.27 2003-03-21 00:52:42 geuzaine Exp $
+// $Id: Simplex.cpp,v 1.28 2003-12-07 00:23:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2003 C. Geuzaine, J.-F. Remacle
 //
@@ -44,6 +44,7 @@ Simplex::Simplex()
   V[0] = V[1] = V[2] = V[3] = NULL;
   S[0] = S[1] = S[2] = S[3] = NULL;
   iEnt = -1;
+  iPart = -1;
   Quality = 0.;
   Num = TotalNumber;
   Visible = VIS_MESH;
@@ -60,6 +61,7 @@ Simplex::Simplex(Vertex * v1, Vertex * v2, Vertex * v3, Vertex * v4)
   Num = TotalNumber;
   THEM->MaxSimplexNum = IMAX(THEM->MaxSimplexNum, Num);
   iEnt = -1;
+  iPart = -1;
   Visible = VIS_MESH;
 }
 
diff --git a/Mesh/Simplex.h b/Mesh/Simplex.h
index f95a3c2bfd32f4d0c12893d6d59b25933f4b46e2..b667ea14bfd25d9bf0e026ef93d0fda2e89e31cd 100644
--- a/Mesh/Simplex.h
+++ b/Mesh/Simplex.h
@@ -30,16 +30,17 @@ typedef struct {
 class Simplex{
 
 public:
-  int     Num;           /* Numero                                       */
-  int     iEnt;          /* Entite geometrique                           */
-  char    Visible;       /* Visualization flag                           */
-  Face    F[4];          /* 4 faces                                      */
-  Vertex  **VSUP;        /* noeuds supplem pour les elts de degre eleves */
-  Vertex  *V[4];         /* 4 noeuds                                     */
-  double  Quality;       /* Qualite du simplexe                          */
-  Coord   Center;        /* centre du CC                                 */
-  double  Radius;        /* Rayon du CC                                  */
-  Simplex *S[4];         /* 4 Voisins                                    */
+  int     Num;           // Number
+  int     iEnt;          // Elementary geometrical entity
+  int     iPart;         // Mesh partition index
+  char    Visible;       // Visualization flag
+  Face    F[4];          // 4 faces
+  Vertex  **VSUP;        // suppl. nodes for higher order elts
+  Vertex  *V[4];         // 4 nodes
+  double  Quality;       // simplex quality
+  Coord   Center;        // CC center
+  double  Radius;        // CC radius
+  Simplex *S[4];         // 4 neighbours
   static  int TotalNumber;
   static  int TotalAllocated;
   Simplex();