diff --git a/DataStr/Tools.cpp b/DataStr/Tools.cpp
index 6359220c8c3c15785ae0af49384716f26e4e75d6..00c06ccb9c4a0de73d783c68bc13f9214b314ece 100644
--- a/DataStr/Tools.cpp
+++ b/DataStr/Tools.cpp
@@ -1,4 +1,4 @@
-// $Id: Tools.cpp,v 1.13 2005-01-01 19:35:27 geuzaine Exp $
+// $Id: Tools.cpp,v 1.14 2005-05-27 19:35:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -23,6 +23,10 @@
 #include <math.h>
 #include "Tools.h"
 
+static List_T *pListeTransfert;
+static Tree_T *pTreeTransfert;
+static Tree_T *pTreeTransfert2;
+
 // Comparison functions
 
 int fcmp_int(const void *a, const void *b)
@@ -65,8 +69,6 @@ List_T *ListOfDouble2ListOfInt(List_T *dList)
 
 // Tree -> List transfer
 
-List_T *pListeTransfert;
-
 void TransfereListe(void *a, void *b)
 {
   List_Add(pListeTransfert, a);
@@ -83,10 +85,21 @@ List_T *Tree2List(Tree_T * pTree)
   return (pListeTransfert);
 }
 
-// Algebraic utilities
+// List -> Tree transfer
+
+void TransfereTree(void *a, void *b)
+{
+  Tree_Add(pTreeTransfert, a);
+}
+
+Tree_T *List2Tree(List_T * pList, int (*fcmp) (const void *a, const void *b))
+{
+  pTreeTransfert = Tree_Create(pList->size, fcmp);
+  List_Action(pList, TransfereTree);
+  return (pTreeTransfert);
+}
 
-Tree_T *pTreeTransfert;
-Tree_T *pTreeTransfert2;
+// Algebraic utilities
 
 void DupliqueArbre(void *a, void *b)
 {
diff --git a/DataStr/Tools.h b/DataStr/Tools.h
index 6838330566bfc3f1cf1659026d064f81cb62713b..155721114e22d19e51a6c6f618f4e6629321e45b 100644
--- a/DataStr/Tools.h
+++ b/DataStr/Tools.h
@@ -29,6 +29,7 @@ int fcmp_double(const void *a, const void *b);
 
 List_T *ListOfDouble2ListOfInt(List_T *dList);
 List_T *Tree2List(Tree_T *pTree);
+Tree_T *List2Tree(List_T * pList, int (*fcmp) (const void *a, const void *b));
 
 Tree_T *Tree_Duplique(Tree_T *pTree);
 Tree_T *Tree_Union(Tree_T *pTreeA, Tree_T *pTreeB);
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 57680a211a6d62c72f669c0bdbe62d45ca7c626f..1f11cfc94191a9148c579596db000e3ac6967127 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.127 2005-05-21 04:55:59 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.128 2005-05-27 19:35:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -77,11 +77,23 @@ void draw_polygon_2d(double r, double g, double b, int n,
 
 int getFirstPhysical(int type, int num)
 {
-  for(int i = 0; i < List_Nbr(THEM->PhysicalGroups); i++){
-    PhysicalGroup *p = *(PhysicalGroup**)List_Pointer(THEM->PhysicalGroups, i);
-    if(p->Typ == type && List_Query(p->Entities, &num, fcmp_int))
-      return p->Num;
-  }  
+  // If we visualize the mesh by coloring the physical entities, this
+  // routine returns the number of the first physical entity of type
+  // "type" that contains the elementary entity "num"
+  if(CTX.mesh.color_carousel == 2){
+    static int warn = 1;
+    if(List_Nbr(THEM->PhysicalGroups) > 100 && warn){
+      Msg(WARNING, "There are many physical entities in the mesh (%d)",
+	  List_Nbr(THEM->PhysicalGroups));
+      Msg(WARNING, "You might want to color the mesh by elementary entity instead");
+      warn = 0;
+    }
+    for(int i = 0; i < List_Nbr(THEM->PhysicalGroups); i++){
+      PhysicalGroup *p = *(PhysicalGroup**)List_Pointer(THEM->PhysicalGroups, i);
+      if(p->Typ == type && List_Query(p->Entities, &num, fcmp_int))
+	return p->Num;
+    }
+  }
   return 0;
 }
 
@@ -299,7 +311,6 @@ void Draw_Mesh_Volume(void *a, void *b)
 
   theVolume = v;
   theColor = v->Color;
-  thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, v->Num);
 
   // we don't use vertex arrays for every volume primitive: only for
   // volume cuts drawn "as surfaces" (using vertex arrays for
@@ -308,6 +319,7 @@ void Draw_Mesh_Volume(void *a, void *b)
      CTX.mesh.vertex_arrays){
     if(CTX.mesh.changed){
       Msg(DEBUG, "regenerate volume mesh vertex arrays");
+      thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, v->Num);
       // triangles
       if(v->TriVertexArray) delete v->TriVertexArray;
       v->TriVertexArray = new VertexArray(3, 1000);
@@ -344,6 +356,7 @@ void Draw_Mesh_Volume(void *a, void *b)
      CTX.mesh.dual || CTX.mesh.volumes_num || CTX.mesh.points_per_element ||
      CTX.mesh.normals){
     Msg(DEBUG, "classic volume data path");
+    thePhysical = getFirstPhysical(MSH_PHYSICAL_VOLUME, v->Num);
     Tree_Action(v->Simplexes, Draw_Mesh_Tetrahedron);
     Tree_Action(v->SimplexesBase, Draw_Mesh_Tetrahedron);
     Tree_Action(v->Hexahedra, Draw_Mesh_Hexahedron);
@@ -362,7 +375,6 @@ void Draw_Mesh_Surface(void *a, void *b)
 
   theSurface = s;
   theColor = s->Color;
-  thePhysical = getFirstPhysical(MSH_PHYSICAL_SURFACE, s->Num);
 
   if(CTX.mesh.changed && CTX.mesh.smooth_normals){
     Msg(DEBUG, "pre-processing smooth normals");
@@ -378,6 +390,7 @@ void Draw_Mesh_Surface(void *a, void *b)
   if(CTX.mesh.vertex_arrays){
     if(CTX.mesh.changed){
       Msg(DEBUG, "regenerate surface mesh vertex arrays");
+      thePhysical = getFirstPhysical(MSH_PHYSICAL_SURFACE, s->Num);
       // triangles
       if(s->TriVertexArray) delete s->TriVertexArray;
       s->TriVertexArray = new VertexArray(3, Tree_Nbr(s->Simplexes) + 
@@ -410,6 +423,7 @@ void Draw_Mesh_Surface(void *a, void *b)
   if(!s->TriVertexArray || CTX.mesh.dual || CTX.mesh.surfaces_num ||
      CTX.mesh.points_per_element || CTX.mesh.normals){
     Msg(DEBUG, "classic triangle data path");
+    thePhysical = getFirstPhysical(MSH_PHYSICAL_SURFACE, s->Num);
     Tree_Action(s->Simplexes, Draw_Mesh_Triangle);
     Tree_Action(s->SimplexesBase, Draw_Mesh_Triangle);
   }
@@ -417,6 +431,7 @@ void Draw_Mesh_Surface(void *a, void *b)
   if(!s->QuadVertexArray || CTX.mesh.dual || CTX.mesh.surfaces_num ||
      CTX.mesh.points_per_element || CTX.mesh.normals){
     Msg(DEBUG, "classic quadrangle data path");
+    thePhysical = getFirstPhysical(MSH_PHYSICAL_SURFACE, s->Num);
     Tree_Action(s->Quadrangles, Draw_Mesh_Quadrangle);
   }
 
@@ -442,11 +457,11 @@ void Draw_Mesh_Curve(void *a, void *b)
 
   theCurve = c;
   theColor = c->Color;
-  thePhysical = getFirstPhysical(MSH_PHYSICAL_LINE, c->Num);
 
   if(CTX.mesh.vertex_arrays){
     if(CTX.mesh.changed){
       Msg(DEBUG, "regenerate curve mesh vertex array");
+      thePhysical = getFirstPhysical(MSH_PHYSICAL_LINE, c->Num);
       if(c->LinVertexArray) delete c->LinVertexArray;
       c->LinVertexArray = new VertexArray(2, Tree_Nbr(c->Simplexes) + 
 					  Tree_Nbr(c->SimplexesBase));
@@ -465,6 +480,7 @@ void Draw_Mesh_Curve(void *a, void *b)
   if(!c->LinVertexArray || CTX.mesh.lines_num ||
      CTX.mesh.points_per_element || CTX.mesh.tangents){
     Msg(DEBUG, "classic line data path");
+    thePhysical = getFirstPhysical(MSH_PHYSICAL_LINE, c->Num);
     Tree_Action(c->Simplexes, Draw_Mesh_Line);
     Tree_Action(c->SimplexesBase, Draw_Mesh_Line);
   }
diff --git a/Mesh/Read_Mesh.cpp b/Mesh/Read_Mesh.cpp
index ed16b303c2199af0c6f4126437694fff1183bacd..f2fac40a02ce52c4e327eb810edb761b2e87219d 100644
--- a/Mesh/Read_Mesh.cpp
+++ b/Mesh/Read_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Read_Mesh.cpp,v 1.88 2005-05-21 01:10:47 geuzaine Exp $
+// $Id: Read_Mesh.cpp,v 1.89 2005-05-27 19:35:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -50,21 +50,6 @@ extern Context_T CTX;
 
 #define NB_NOD_MAX_ELM 30
 
-void addPhysicalGroup(Mesh * M, int Type, int Physical, int Elementary)
-{
-  PhysicalGroup *pg;
-  if((pg = FindPhysicalGroup(Physical, Type, M))) {
-    List_Insert(pg->Entities, &Elementary, fcmp_int);
-  }
-  else {
-    List_T *tmp = List_Create(1, 1, sizeof(int));
-    List_Add(tmp, &Elementary);
-    pg = Create_PhysicalGroup(Physical, Type, tmp);
-    List_Add(M->PhysicalGroups, &pg);
-    List_Delete(tmp);
-  }
-}
-
 // If a "normal" elementary entity does not exist, we create a
 // "discrete" entity, i.e., an entity entirely defined by its mesh
 
@@ -98,6 +83,45 @@ Volume *addElementaryVolume(Mesh * M, int Num)
   return v;
 }
 
+void addPhysicalGroup(Tree_T * groups, int Type, int Physical, int Elementary)
+{
+  // we add in a temporary group tree for performance reasons; the
+  // tree is converted back into a list in the mesh at the end of read_mesh
+  PhysicalGroup g, *pg;
+  pg = &g;
+  pg->Num = Physical;
+  pg->Typ = Type;
+  if(Tree_Query(groups, &pg)) {
+    List_Insert(pg->Entities, &Elementary, fcmp_int);
+  }
+  else {
+    List_T *tmp = List_Create(1, 1, sizeof(int));
+    List_Add(tmp, &Elementary);
+    pg = Create_PhysicalGroup(Physical, Type, tmp);
+    Tree_Add(groups, &pg);
+    List_Delete(tmp);
+  }
+}
+
+int addMeshPartition(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;
+  }
+  return 0;
+}
+
 int getNbrNodes(int Type)
 {
   switch (Type) {
@@ -154,6 +178,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
   Surface *s;
   Volume *v;
   Tree_T *Duplicates = NULL;
+  Tree_T *groups = List2Tree(M->PhysicalGroups, comparePhysicalGroup);
 
   while(1) {
     do {
@@ -262,7 +287,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 	if(version <= 1.0){
 	  fscanf(fp, "%d %d %d %d %d",
 		 &Num, &Type, &Physical, &Elementary, &Nbr_Nodes);
-	  Partition = Physical;
+	  Partition = 1;
 	  int Nbr_Nodes_Check = getNbrNodes(Type);
 	  if(!Nbr_Nodes_Check){
 	    Msg(GERROR, "Unknown type for element %d", Num); 
@@ -325,11 +350,11 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case LGN1:
         case LGN2:
 	  c = addElementaryCurve(M, abs(Elementary));
-	  addPhysicalGroup(M, MSH_PHYSICAL_LINE, Physical, abs(Elementary));
+	  addPhysicalGroup(groups, MSH_PHYSICAL_LINE, Physical, abs(Elementary));
           simp = Create_SimplexBase(vertsp[0], vertsp[1], NULL, NULL);
           simp->Num = Num;
           simp->iEnt = Elementary;
-          simp->iPart = Add_MeshPartition(Partition, M);
+          simp->iPart = addMeshPartition(Partition, M);
 	  if(Type == LGN2){
 	    simp->VSUP = (Vertex **) Malloc(1 * sizeof(Vertex *));
 	    simp->VSUP[0] = vertsp[2];
@@ -355,11 +380,11 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case TRI1:
         case TRI2:
 	  s = addElementarySurface(M, Elementary);
-	  addPhysicalGroup(M, MSH_PHYSICAL_SURFACE, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_SURFACE, Physical, Elementary);
           simp = Create_SimplexBase(vertsp[0], vertsp[1], vertsp[2], NULL);
           simp->Num = Num;
           simp->iEnt = Elementary;
-          simp->iPart = Add_MeshPartition(Partition, M);
+          simp->iPart = addMeshPartition(Partition, M);
 	  if(Type == TRI2){
 	    simp->VSUP = (Vertex **) Malloc(3 * sizeof(Vertex *));
 	    for(i = 0; i < 3; i++){
@@ -379,11 +404,11 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case QUA1:
         case QUA2:
 	  s = addElementarySurface(M, Elementary);
-	  addPhysicalGroup(M, MSH_PHYSICAL_SURFACE, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_SURFACE, Physical, Elementary);
           quad = Create_Quadrangle(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
           quad->Num = Num;
           quad->iEnt = Elementary;
-          quad->iPart = Add_MeshPartition(Partition, M);
+          quad->iPart = addMeshPartition(Partition, M);
 	  if(Type == QUA2){
 	    quad->VSUP = (Vertex **) Malloc((4 + 1) * sizeof(Vertex *));
 	    for(i = 0; i < 4 + 1; i++){
@@ -403,11 +428,11 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case TET1:
         case TET2:
 	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(M, MSH_PHYSICAL_VOLUME, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
           simp = Create_SimplexBase(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
           simp->Num = Num;
           simp->iEnt = Elementary;
-          simp->iPart = Add_MeshPartition(Partition, M);
+          simp->iPart = addMeshPartition(Partition, M);
 	  if(Type == TET2){
 	    simp->VSUP = (Vertex **) Malloc(6 * sizeof(Vertex *));
 	    for(i = 0; i < 6; i++){
@@ -429,12 +454,12 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case HEX1:
         case HEX2:
 	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(M, MSH_PHYSICAL_VOLUME, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
           hex = Create_Hexahedron(vertsp[0], vertsp[1], vertsp[2], vertsp[3],
                                   vertsp[4], vertsp[5], vertsp[6], vertsp[7]);
           hex->Num = Num;
           hex->iEnt = Elementary;
-          hex->iPart = Add_MeshPartition(Partition, M);
+          hex->iPart = addMeshPartition(Partition, M);
 	  if(Type == HEX2){
 	    hex->VSUP = (Vertex **) Malloc((12 + 6 + 1) * sizeof(Vertex *));
 	    for(i = 0; i < 12 + 6 + 1; i++){
@@ -456,12 +481,12 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case PRI1:
         case PRI2:
 	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(M, MSH_PHYSICAL_VOLUME, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
           pri = Create_Prism(vertsp[0], vertsp[1], vertsp[2],
                              vertsp[3], vertsp[4], vertsp[5]);
           pri->Num = Num;
           pri->iEnt = Elementary;
-          pri->iPart = Add_MeshPartition(Partition, M);
+          pri->iPart = addMeshPartition(Partition, M);
 	  if(Type == PRI2){
 	    pri->VSUP = (Vertex **) Malloc((9 + 3) * sizeof(Vertex *));
 	    for(i = 0; i < 9 + 3; i++){
@@ -483,12 +508,12 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
         case PYR1:
         case PYR2:
 	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(M, MSH_PHYSICAL_VOLUME, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
           pyr = Create_Pyramid(vertsp[0], vertsp[1], vertsp[2],
                                vertsp[3], vertsp[4]);
           pyr->Num = Num;
           pyr->iEnt = Elementary;
-          pyr->iPart = Add_MeshPartition(Partition, M);
+          pyr->iPart = addMeshPartition(Partition, M);
 	  if(Type == PYR2){
 	    pyr->VSUP = (Vertex **) Malloc((8 + 1) * sizeof(Vertex *));
 	    for(i = 0; i < 8 + 1; i++){
@@ -508,7 +533,7 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 #endif
           break;
         case PNT:
-	  addPhysicalGroup(M, MSH_PHYSICAL_POINT, Physical, Elementary);
+	  addPhysicalGroup(groups, MSH_PHYSICAL_POINT, Physical, Elementary);
 	  // we need to make a new one: vertices in M->Vertices and
 	  // M->Points should never point to the same memory location
 	  vert = Create_Vertex(vertsp[0]->Num, vertsp[0]->Pos.X, vertsp[0]->Pos.Y, 
@@ -562,6 +587,11 @@ void Read_Mesh_MSH(Mesh * M, FILE * fp)
 
   // Transfer the vertices from temp trees into lists
   Tree_Action(M->Curves, Transfer_VertexTree2List);
+
+  // Transfer the temp group tree back into the mesh
+  List_Delete(M->PhysicalGroups);
+  M->PhysicalGroups = Tree2List(groups);
+  Tree_Delete(groups);
 }
 
 // Read mesh in VTK format