Index: gmsh_io.C
===================================================================
--- gmsh_io.C	(revision 3513)
+++ gmsh_io.C	(working copy)
@@ -339,6 +339,7 @@
   char       buf[bufLen+1];
   int        format=0, size=0;
   Real       version = 1.0;
+  std::vector<subdomain_id_type> sbd_ids ;
   
   // map to hold the node numbers for translation
   // note the the nodes can be non-consecutive
@@ -412,7 +413,7 @@
           for (unsigned int iel=0; iel<numElem; ++iel)      
             {
               unsigned int id, type, physical, elementary,
-                /* partition = 1,*/ nnodes, ntags;
+                partition = 1, nnodes, ntags;
 	      // note - partition was assigned but never used - BSK
               if(version <= 1.0)
                 {
@@ -430,12 +431,12 @@
                         physical = tag;
                       else if(j == 1)
                         elementary = tag;
-                      // else if(j == 2)
-                      //  partition = tag;
+                      else if(j == 2)
+                       partition = tag;
                       // ignore any other tags for now
                     }
                 }
-
+				
               // consult the import element table which element to build
               const elementDefinition& eletype = eletypes_imp[type];
               nnodes = eletype.nnodes;
@@ -448,6 +449,22 @@
                   // add the elements to the mesh
                   Elem* elem = Elem::build(eletype.type).release();
                   elem->set_id(elem_id_counter);
+				  
+				  subdomain_id_type& sbdid = elem->subdomain_id() ;
+				  sbdid = static_cast<subdomain_id_type>(physical) ;
+                  
+				  // we can save the elementary surface/volume ids
+				  // for each element, if physics can use this information.
+				  // volume_id_type& volid = elem->volume_id() ;
+				  // volid = static_cast<volume_id_type>(elementary) ;
+				  
+				  // search if the vector has this subdomain id added before; if not add it
+				  std::vector<subdomain_id_type>::const_iterator sbd_it = std::find(sbd_ids.begin(), sbd_ids.end(), sbdid) ;
+				  if (sbd_it == sbd_ids.end())
+				  {
+					sbd_ids.push_back(sbdid);
+				  }
+					
                   mesh.add_elem(elem);
 
 		  // different to iel, lower dimensional elems aren't added
@@ -487,7 +504,7 @@
                     }
 
                     // Finally, set the subdomain ID to physical
-                    elem->subdomain_id() = physical;
+                    elem->subdomain_id() = sbdid;
                 } // if element.dim == dim
               // if this is a boundary
               else if (eletype.dim == dim-1)
@@ -586,6 +603,9 @@
     } // while !in.eof()
 
   }
+  
+  unsigned int& nsubdms = mesh.set_n_subdomains() ;
+  nsubdms = sbd_ids.size() ;
 
 }
 
@@ -681,8 +701,14 @@
         out << eletype.exptype;
 
         // write the number of tags and
-        // tag1 (physical entity), and tag2 (geometric entity)
-        out << " 3 1 1 ";
+        // tag1 (physical entity), and tag2 (geometric entity), and tag3 (partition entity)
+        out << " 3 ";
+		
+		out << elem->subdomain_id() << " " ;
+		
+		// this should be surface/volume id, if we had saved this 
+		// for each elem while reading the mesh
+		out << elem->id() << " " ;
 
         // write the partition the element belongs to
         out << elem->processor_id()+1 << " ";