diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp
index 212dbc6c..8e656938 100644
--- a/libsrc/meshing/parallelmesh.cpp
+++ b/libsrc/meshing/parallelmesh.cpp
@@ -31,10 +31,9 @@ namespace netgen
     if (id == 0)
       PrintMessage (1, "Send/Receive mesh");
 
-
     {
       // distribute global information
-      int nelglob, nvglob;
+      int nelglob, nseglob, nvglob;
       if (id == 0)
 	{
 	  paralleltop -> SetNV (GetNV());
@@ -43,15 +42,19 @@ namespace netgen
 	  paralleltop -> SetNSE (GetNSE());
 	  
 	  nelglob = GetNE();
+	  nseglob = GetNSE();
 	  nvglob = GetNV();
 	}
       
+      MyMPI_Bcast (dimension);
       MyMPI_Bcast (nelglob);
+      MyMPI_Bcast (nseglob);
       MyMPI_Bcast (nvglob);
       
       if (id > 0)
 	{
 	  paralleltop -> SetNEGlob (nelglob);
+	  paralleltop -> SetNSEGlob (nseglob);
 	  paralleltop -> SetNVGlob (nvglob);
 	}
     }
@@ -98,25 +101,31 @@ namespace netgen
   {
     Array<MPI_Request> sendrequests;
 
+    PrintMessage ( 3, "Sending vertices");
+
+
     Array<int> num_els_on_proc(ntasks);
     num_els_on_proc = 0;
     for (ElementIndex ei = 0; ei < GetNE(); ei++)
       num_els_on_proc[(*this)[ei].GetPartition()]++;
 
-
-    // get number of vertices for each processor
-    Array<int> nelloc (ntasks);
-    
-    nelloc = 0;
-    // elarraysize = 1;
-
-    PrintMessage ( 3, "Sending vertices");
-
     TABLE<ElementIndex> els_of_proc (num_els_on_proc);
     for (ElementIndex ei = 0; ei < GetNE(); ei++)
       els_of_proc.Add ( (*this)[ei].GetPartition(), ei);
 
 
+    Array<int> num_sels_on_proc(ntasks);
+    num_sels_on_proc = 0;
+    for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
+      num_sels_on_proc[(*this)[ei].GetPartition()]++;
+
+    TABLE<SurfaceElementIndex> sels_of_proc (num_sels_on_proc);
+    for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
+      sels_of_proc.Add ( (*this)[ei].GetPartition(), ei);
+
+
+
+
     Array<int, PointIndex::BASE> vert_flag (GetNV());
     Array<int, PointIndex::BASE> num_procs_on_vert (GetNV());
     Array<int> num_verts_on_proc (ntasks);
@@ -125,6 +134,12 @@ namespace netgen
     num_procs_on_vert = 0;
     vert_flag = -1;
 
+    Array<int> nelloc (ntasks);
+    nelloc = 0;
+    Array<int> nselloc (ntasks);
+    nselloc = 0;
+
+
     for (int dest = 1; dest < ntasks; dest++)
       {
 	FlatArray<ElementIndex> els = els_of_proc[dest];
@@ -145,11 +160,32 @@ namespace netgen
 		    paralleltop -> SetDistantPNum ( dest, epi, num_verts_on_proc[dest]);
 		  }
 	      }
-
-	    // elarraysize[dest] += 3 + el.GetNP();
 	    nelloc[dest] ++;
 	    paralleltop -> SetDistantEl ( dest, els[hi]+1, nelloc[dest] );
 	  }
+
+
+	FlatArray<SurfaceElementIndex> sels = sels_of_proc[dest];
+
+	for (int hi = 0; hi < sels.Size(); hi++)
+	  {
+	    const Element2d & el = (*this) [ sels[hi] ];
+	    for (int i = 0; i < el.GetNP(); i++)
+	      {
+		PointIndex epi = el[i]; 
+		if (vert_flag[epi] < dest)
+		  {
+		    vert_flag[epi] = dest;
+
+		    num_verts_on_proc[dest]++;
+		    num_procs_on_vert[epi]++;
+
+		    paralleltop -> SetDistantPNum ( dest, epi, num_verts_on_proc[dest]);
+		  }
+	      }
+	    nselloc[dest] ++;
+	    paralleltop -> SetDistantSurfEl ( dest, sels[hi]+1, nselloc[dest] );
+	  }
       }
 
     TABLE<PointIndex> verts_of_proc (num_verts_on_proc);
@@ -163,7 +199,6 @@ namespace netgen
 	for (int hi = 0; hi < els.Size(); hi++)
 	  {
 	    const Element & el = (*this) [ els[hi] ];
-	    
 	    for (int i = 0; i < el.GetNP(); i++)
 	      {
 		PointIndex epi = el[i];
@@ -174,6 +209,22 @@ namespace netgen
 		  }
 	      }
 	  }
+
+	FlatArray<SurfaceElementIndex> sels = sels_of_proc[dest];
+	for (int hi = 0; hi < sels.Size(); hi++)
+	  {
+	    const Element2d & el = (*this) [ sels[hi] ];
+	    for (int i = 0; i < el.GetNP(); i++)
+	      {
+		PointIndex epi = el[i];
+		if (vert_flag[epi] < dest)
+		  {
+		    vert_flag[epi] = dest;
+		    procs_of_vert.Add (epi, dest);
+		  }
+	      }
+	  }
+
       }
 
     for (int vert = 1; vert <= GetNP(); vert++ )
@@ -241,6 +292,8 @@ namespace netgen
     for ( int dest = 1; dest < ntasks; dest ++ )
       sendrequests.Append (MyMPI_ISend (distpnums[dest], dest, MPI_TAG_MESH+1));
 
+
+
     PrintMessage ( 3, "Sending elements" );
 
     Array<int> elarraysize (ntasks);
@@ -270,7 +323,6 @@ namespace netgen
       sendrequests.Append (MyMPI_ISend (elementarrays[dest], dest, MPI_TAG_MESH+2));
 
     PrintMessage ( 3, "Sending Face Descriptors" );
-  
 
     Array<double> fddata (6 * GetNFD());
     for (int fdi = 1; fdi <= GetNFD(); fdi++)
@@ -288,71 +340,37 @@ namespace netgen
 
     PrintMessage ( 3, "Sending Surface elements" );
   
-    Array <int> nlocsel(ntasks), bufsize(ntasks); // , seli(ntasks);
-    for ( int i = 0; i < ntasks; i++)
-      {
-	nlocsel[i] = 0;
-	bufsize[i] = 1;
-      }
+    Array <int> nlocsel(ntasks), bufsize(ntasks); 
+    nlocsel = 0;
+    bufsize = 1;
 
     for (int sei = 1; sei <= GetNSE(); sei++ )
       {
-	int ei1, ei2;
-	GetTopology().GetSurface2VolumeElement (sei, ei1, ei2);
 	const Element2d & sel = SurfaceElement (sei);
-
-        for (int j = 0; j < 2; j++)
-          {
-            int ei = (j == 0) ? ei1 : ei2;
-            if ( ei > 0 && ei <= GetNE() )
-              {
-                const Element & el = VolumeElement (ei);
-                int dest = el.GetPartition();
-		nlocsel[dest] ++;
-		bufsize[dest] += 4 + 2*sel.GetNP();
-              }
-	  }
+	int dest = sel.GetPartition();
+	nlocsel[dest] ++;
+	bufsize[dest] += 4 + 2*sel.GetNP();
       }
     
     TABLE<int> selbuf(bufsize);
 
-    Array<int> nselloc (ntasks);
-    nselloc = 0;
-
     for (int dest = 1; dest < ntasks; dest++ )
       selbuf.Add (dest, nlocsel[dest]);
 
     for (int sei = 1; sei <= GetNSE(); sei ++ )
       {
-	int ei1, ei2;
-	GetTopology().GetSurface2VolumeElement (sei, ei1, ei2);
 	const Element2d & sel = SurfaceElement (sei);
+	int dest = sel.GetPartition();
 
-	int isghost = 0;
-        for (int j = 0; j < 2; j++)
-          {
-            int ei = (j == 0) ? ei1 : ei2;
-            if ( ei > 0 && ei <= GetNE() )
-              {
-                const Element & el = VolumeElement (ei);
-                int dest = el.GetPartition();
-                if (dest > 0)
-                  {
-		    selbuf.Add (dest, sei);
-		    selbuf.Add (dest, sel.GetIndex());
-		    selbuf.Add (dest, isghost);
-		    selbuf.Add (dest, sel.GetNP());
-
-                    for ( int ii = 1; ii <= sel.GetNP(); ii++)
-                      {
-                        selbuf.Add (dest, sel.PNum(ii));
-                        selbuf.Add (dest, sel.GeomInfoPi(ii).trignum);
-                      }
-		    nselloc[dest] ++;
-		    paralleltop -> SetDistantSurfEl ( dest, sei, nselloc[dest] );
-		    isghost = 1;
-                  }
-              }
+	selbuf.Add (dest, sei);
+	selbuf.Add (dest, sel.GetIndex());
+	selbuf.Add (dest, 0);
+	selbuf.Add (dest, sel.GetNP());
+	
+	for ( int ii = 1; ii <= sel.GetNP(); ii++)
+	  {
+	    selbuf.Add (dest, sel.PNum(ii));
+	    selbuf.Add (dest, sel.GeomInfoPi(ii).trignum);
 	  }
       }
 
@@ -373,13 +391,13 @@ namespace netgen
       {
 	Array<int> volels;
 	const MeshTopology & topol = GetTopology();
-	topol . GetSegmentVolumeElements ( segi, volels );
+	topol . GetSegmentSurfaceElements ( segi, volels );
         for (int j = 0; j < volels.Size(); j++)
           {
             int ei = volels[j];
-            if ( ei > 0 && ei <= GetNE() )
+            if ( ei > 0 && ei <= GetNSE() )
               {
-                const Element & el = VolumeElement (ei);
+                const Element2d & el = SurfaceElement (ei);
                 int dest = el.GetPartition();
 		nlocseg[dest] ++;
 		bufsize[dest] += 14;
@@ -392,15 +410,15 @@ namespace netgen
     for ( int ls=1; ls <= GetNSeg(); ls++)
       {
 	Array<int> volels;
-	GetTopology().GetSegmentVolumeElements ( ls, volels );
+	GetTopology().GetSegmentSurfaceElements ( ls, volels );
 	const Segment & seg = LineSegment (ls);
 
         for (int j = 0; j < volels.Size(); j++)
           {
             int ei = volels[j];
-            if ( ei > 0 && ei <= GetNE() )
+            if ( ei > 0 && ei <= GetNSE() )
               {
-                const Element & el = VolumeElement (ei);
+                const Element2d & el = SurfaceElement (ei);
                 int dest = el.GetPartition();
 
 		if ( dest > 0 )
@@ -430,6 +448,8 @@ namespace netgen
       sendrequests.Append (MyMPI_ISend(segmbuf[dest], dest, MPI_TAG_MESH+5));
 
     MPI_Waitall (sendrequests.Size(), &sendrequests[0], MPI_STATUS_IGNORE);
+
+    MPI_Barrier(MPI_COMM_WORLD);
   }
 
 
@@ -451,8 +471,6 @@ namespace netgen
 
     // string st;
 
-
-    
     // receive vertices
     NgProfiler::StartTimer (timer_pts);
     
@@ -463,7 +481,7 @@ namespace netgen
     paralleltop -> SetNV (numvert);
     
     // INDEX_CLOSED_HASHTABLE<int> glob2loc_vert_ht (3*numvert+1);
-    INDEX_HASHTABLE<int> glob2loc_vert_ht (100);
+    INDEX_HASHTABLE<int> glob2loc_vert_ht (3*numvert+1);
 
     for (int vert = 0; vert < numvert; vert++)
       {
@@ -488,7 +506,6 @@ namespace netgen
     
     NgProfiler::StopTimer (timer_pts);
     *testout << "got " << numvert << " vertices" << endl;
-    
 
     {
       Element el;
@@ -498,7 +515,6 @@ namespace netgen
       
       NgProfiler::RegionTimer reg(timer_els);
 
-      // double sumcosts = 0, sum = 0;
       for (int ind = 0, elnum = 1; ind < elarray.Size(); elnum++)
 	{
 	  paralleltop->SetLoc2Glob_VolEl ( elnum,  elarray[ind++]);
@@ -506,27 +522,13 @@ namespace netgen
 	  el.SetIndex(elarray[ind++]);
 	  el.SetNP(elarray[ind++]);
 	  
-	  /*
-	  for ( int j = 0; j < el.GetNP(); j++)
-	    {
-	      int costs = 1; // glob2loc_vert_ht.CalcPositionCosts (elarray[ind+j]); 
-	      // *testout << "costs to hash " << elarray[ind+j] << ", hashvalue = " 
-	      // << glob2loc_vert_ht.HashValue(elarray[ind+j]) << ": " << costs << endl;
-	      sumcosts += costs;
-	      sum++;
-	    }
-	  */
-
 	  for ( int j = 0; j < el.GetNP(); j++)
 	    el[j] = glob2loc_vert_ht.Get (elarray[ind++]); 
 	  
 	  AddVolumeElement (el);
 	}
-      // *testout << "average hashing costs: " << sumcosts / sum << endl;
-      // cout << "average hashing costs: " << sumcosts / sum << endl;
     }
 
-
     {
       Array<double> fddata;
       MyMPI_Recv (fddata, 0, MPI_TAG_MESH+3);
@@ -541,7 +543,6 @@ namespace netgen
 
     {
       NgProfiler::RegionTimer reg(timer_sels);
-      // surface elements
       Array<int> selbuf;
 
       MyMPI_Recv ( selbuf, 0, MPI_TAG_MESH+4);
@@ -620,12 +621,14 @@ namespace netgen
 	}
     }
     
+    MPI_Barrier(MPI_COMM_WORLD);
 
     int timerloc = NgProfiler::CreateTimer ("Update local mesh");
     int timerloc2 = NgProfiler::CreateTimer ("CalcSurfacesOfNode");
 
     NgProfiler::RegionTimer regloc(timerloc);
     PrintMessage (2, "Got ", GetNE(), " elements");
+    PrintMessage (2, "Got ", GetNSE(), " surface elements");
 
     NgProfiler::StartTimer (timerloc2);
 
@@ -638,6 +641,7 @@ namespace netgen
     
     SetNextMajorTimeStamp();
     // paralleltop->Print();
+    cout << "receive mesh complete, id = " << id << endl;
   }
   
 
@@ -648,14 +652,10 @@ namespace netgen
   // call it only for the master !
   void Mesh :: Distribute ()
   {
-    if ( id != 0 || ntasks == 1 ) return;
-    // metis partition of mesh, only if more than one proc
-    
+    if (id != 0 || ntasks == 1 ) return;
 
 
-    // partition mesh
 #ifdef METIS
-
     ParallelMetis ();
 #else
     for (ElementIndex ei = 0; ei < GetNE(); ei++)
@@ -668,17 +668,11 @@ namespace netgen
     for (SurfaceElementIndex ei = 0; ei < GetNSE(); ei++)
       *testout << "sel(" << int(ei) << ") is in part " << (*this)[ei].GetPartition() << endl;
 
-
     // send partition
-
-    //    for (int dest = 1; dest < ntasks; dest++)
-    //      MyMPI_Send ("mesh", dest, MPI_TAG_CMD);
     MyMPI_SendCmd ("mesh");
     SendRecvMesh (); 
 
-
     paralleltop -> UpdateCoarseGrid();
-
     // paralleltop -> Print();
   }
   
@@ -704,50 +698,42 @@ namespace netgen
     int ne = GetNE();
     int nn = GetNP();
 
-    if (ntasks <= 2)
+    if (ntasks <= 2 || ne <= 1)
       {
         if (ntasks == 1) return;
         
         for (int i=1; i<=ne; i++)
           VolumeElement(i).SetPartition(1);
 
+        for (int i=1; i<=GetNSE(); i++)
+          SurfaceElement(i).SetPartition(1);
+
         return;
       }
 
 
     bool uniform_els = true;
-    ELEMENT_TYPE elementtype = TET; // VolumeElement(1).GetType();
-    // metis works only for uniform tet/hex meshes
-    for ( int el = 2; el <= GetNE(); el++ )
+
+    ELEMENT_TYPE elementtype = TET; 
+    for ( int el = 1; el <= GetNE(); el++ )
       if ( VolumeElement(el).GetType() != elementtype )
 	{
-// 	  int nelperproc = ne / (ntasks-1);
-// 	  for (int i=1; i<=ne; i++)
-// 	    { 
-// 	      int partition = i / nelperproc + 1;
-// 	      if ( partition >= ntasks ) partition = ntasks-1;
-// 	      VolumeElement(i).SetPartition(partition);
-// 	    }
-	  
 	  uniform_els = false;
 	  break;
 	}
 
-    // uniform_els = false;
+
     if (!uniform_els)
       {
-	PartHybridMesh ( );   // neloc );
+	PartHybridMesh ( );  
 	return;
       }
 
 
-
     // uniform (TET) mesh,  JS
     int npe = VolumeElement(1).GetNP();
-    
-    idxtype *elmnts;
-    elmnts = new idxtype[ne*npe];
-    
+    Array<idxtype> elmnts(ne*npe);
+
     int etype;
     if (elementtype == TET)
       etype = 2;
@@ -763,10 +749,7 @@ namespace netgen
     int nparts = ntasks-1;
     
     int edgecut;
-    idxtype *epart, *npart;
-    epart = new idxtype[ne];
-    npart = new idxtype[nn];
-    
+    Array<idxtype> epart(ne), npart(nn);
 
 //     if ( ntasks == 1 ) 
 //       {
@@ -787,35 +770,49 @@ namespace netgen
 //       }
     
 
+
     cout << "call metis ... " << flush;
 
     int timermetis = NgProfiler::CreateTimer ("Metis itself");
     NgProfiler::StartTimer (timermetis);
 
-    METIS_PartMeshDual (&ne, &nn, elmnts, &etype, &numflag, &nparts,
-			&edgecut, epart, npart);
+    METIS_PartMeshDual (&ne, &nn, &elmnts[0], &etype, &numflag, &nparts,
+			&edgecut, &epart[0], &npart[0]);
 
     NgProfiler::StopTimer (timermetis);
 
     cout << "complete" << endl;
-    cout << "edge-cut: " << edgecut << ", balance: " << ComputeElementBalance(ne, nparts, epart) << endl;
-    
+    cout << "edge-cut: " << edgecut << ", balance: " 
+	 << ComputeElementBalance(ne, nparts, &epart[0]) << endl;
 
     // partition numbering by metis : 0 ...  ntasks - 1
     // we want:                       1 ...  ntasks
-
     for (int i=1; i<=ne; i++)
       VolumeElement(i).SetPartition(epart[i-1] + 1);
 
 
-    delete [] elmnts; 
-    delete [] epart; 
-    delete [] npart;
+    for (int sei = 1; sei <= GetNSE(); sei++ )
+      {
+	int ei1, ei2;
+	GetTopology().GetSurface2VolumeElement (sei, ei1, ei2);
+	Element2d & sel = SurfaceElement (sei);
+
+        for (int j = 0; j < 2; j++)
+          {
+            int ei = (j == 0) ? ei1 : ei2;
+            if ( ei > 0 && ei <= GetNE() )
+              {
+		sel.SetPartition (VolumeElement(ei).GetPartition());
+		break;
+	      }
+	  }	
+      }
+    
   }
 #endif
 
 
-  void Mesh :: PartHybridMesh () //  Array<int> & neloc ) 
+  void Mesh :: PartHybridMesh () 
   {
 #ifdef METIS
     int ne = GetNE();
diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp
index 73722935..4a03a377 100644
--- a/libsrc/meshing/paralleltop.cpp
+++ b/libsrc/meshing/paralleltop.cpp
@@ -429,16 +429,12 @@ namespace netgen
     if (id == 0)
       PrintMessage ( 3, "UPDATE GLOBAL COARSEGRID STARTS" );      // JS
 
+
     MPI_Group MPI_GROUP_WORLD;
-    
-    int n_ho = netgen::ntasks - 1;
-    int * process_ranks = new int[netgen::ntasks-1];
-    for ( int i = 0; i < netgen::ntasks-1; i++ )
-      process_ranks[i] = i+1;
-    
-    MPI_Comm_group ( MPI_COMM_WORLD, &MPI_GROUP_WORLD);
-    MPI_Group_incl ( MPI_GROUP_WORLD, n_ho, process_ranks, & MPI_HIGHORDER_WORLD);
-    MPI_Comm_create ( MPI_COMM_WORLD, MPI_HIGHORDER_WORLD, & MPI_HIGHORDER_COMM);
+    int process_ranks[] = { 0 };
+    MPI_Comm_group (MPI_COMM_WORLD, &MPI_GROUP_WORLD);
+    MPI_Group_excl (MPI_GROUP_WORLD, 1, process_ranks, &MPI_HIGHORDER_WORLD);
+    MPI_Comm_create (MPI_COMM_WORLD, MPI_HIGHORDER_WORLD, &MPI_HIGHORDER_COMM);
 
 
     int timer = NgProfiler::CreateTimer ("UpdateCoarseGridGlobal");
@@ -446,6 +442,8 @@ namespace netgen
 
 
     *testout << "ParallelMeshTopology :: UpdateCoarseGridGlobal" << endl;
+    cout << "update global, id = " << id << endl;
+
     const MeshTopology & topology = mesh.GetTopology();
   
     Array<int> sendarray, recvarray;
@@ -479,152 +477,82 @@ namespace netgen
 	for ( int i = 0; i < topology .GetNFaces(); i++)
 	  loc2distface[i][0] = i+1;
       }
-  
-
-    if ( id == 0 )
-      sendarray.Append (nfa);
-
-    BitArray recvface(nfa);
-    recvface.Clear();
-   
-    /*
-      Array<int> edges, pnums, faces;
-      for ( int el = 1; el <= ne; el++ )
-      {
-      topology.GetElementFaces (el, faces);
-      int globeli = GetLoc2Glob_VolEl(el);
-
-      for ( int fai = 0; fai < faces.Size(); fai++)
-      {
-      int fa = faces[fai];
-
-      topology.GetElementEdges ( el, edges );
-      topology.GetFaceVertices ( fa, pnums );
-            
-      // send : 
-      // localfacenum
-      // np
-      // ned
-      // globalpnums
-      // localpnums
-            
-      // localedgenums mit globalv1, globalv2
-                
-      sendarray. Append ( fa );
-      sendarray. Append ( globeli );
-      sendarray. Append ( pnums.Size() );
-      sendarray. Append ( edges.Size() );
-
-      if (id == 0)
-      for ( int i = 0; i < pnums.Size(); i++ )
-      sendarray. Append( pnums[i] );
-      else
-      for ( int i = 0; i < pnums.Size(); i++ )
-      sendarray. Append( GetLoc2Glob_Vert(pnums[i]) );
-
-      for ( int i = 0; i < pnums.Size(); i++ )
-      sendarray. Append(pnums[i] );
-                
-      for ( int i = 0; i < edges.Size(); i++ )
-      {
-      sendarray. Append(edges[i] );
-      int v1, v2;
-      topology . GetEdgeVertices ( edges[i], v1, v2 );
-      int dv1 = GetLoc2Glob_Vert ( v1 );
-      int dv2 = GetLoc2Glob_Vert ( v2 );
-                
-      if (id > 0) if ( dv1 > dv2 ) swap ( dv1, dv2 );
-      sendarray . Append ( dv1 );
-      sendarray . Append ( dv2 );
-      }
-      }
-      }
-    */
-
-    // new version
-    Array<int> edges, pnums, faces, elpnums;
-    sendarray.Append (ne);
-    for ( int el = 1; el <= ne; el++ )
-      {
-	topology.GetElementFaces (el, faces);
-	topology.GetElementEdges ( el, edges );
-	const Element & volel = mesh.VolumeElement (el);
-
-	int globeli = GetLoc2Glob_VolEl(el);
-
-	sendarray. Append ( globeli );
-	sendarray. Append ( faces.Size() );
-	sendarray. Append ( edges.Size() );
-	sendarray. Append ( volel.GetNP() );
-
-	for ( int i = 0; i < faces.Size(); i++ )
-	  sendarray. Append(faces[i] );
-	for ( int i = 0; i < edges.Size(); i++ )
-	  sendarray. Append(edges[i] );
-
-	for ( int i = 0; i < volel.GetNP(); i++ )
-	  if (id == 0)
-	    sendarray. Append(volel[i] );
-	  else
-	    sendarray. Append(GetLoc2Glob_Vert (volel[i]));
-      }
-    // end new version
 
 
-
-    BitArray edgeisinit(ned), vertisinit(np);
-    edgeisinit.Clear();
-    vertisinit.Clear();
     
-    // Array for temporary use, to find local from global element fast
-    Array<int,1> glob2loc_el;
-    if ( id != 0 ) 
+    if ( id == 0 )
       {
+	sendarray.Append (nfa);
+	sendarray.Append (ned);
+
+	Array<int> edges, pnums, faces, elpnums;
+	sendarray.Append (ne);
+	for ( int el = 1; el <= ne; el++ )
+	  {
+	    topology.GetElementFaces (el, faces);
+	    topology.GetElementEdges ( el, edges );
+	    const Element & volel = mesh.VolumeElement (el);
+
+	    int globeli = GetLoc2Glob_VolEl(el);
+	    // cout << "el = " << el << ", globeli = " << globeli << endl;
+
+	    sendarray. Append ( globeli );
+	    sendarray. Append ( faces.Size() );
+	    sendarray. Append ( edges.Size() );
+	    sendarray. Append ( volel.GetNP() );
+
+	    for ( int i = 0; i < faces.Size(); i++ )
+	      sendarray. Append(faces[i] );
+	    for ( int i = 0; i < edges.Size(); i++ )
+	      sendarray. Append(edges[i] );
+
+	    for ( int i = 0; i < volel.GetNP(); i++ )
+	      if (id == 0)
+		sendarray. Append(volel[i] );
+	      else
+		sendarray. Append(GetLoc2Glob_Vert (volel[i]));
+	  }
+
+
+	sendarray.Append (nsurfel);
+	for (int el = 1; el <= nsurfel; el++)
+	  {
+	    topology.GetSurfaceElementEdges ( el, edges );
+	    const Element2d & surfel = mesh.SurfaceElement (el);
+
+	    sendarray. Append ( el );
+	    sendarray. Append ( edges.Size() );
+	    sendarray. Append ( surfel.GetNP() );
+
+	    for ( int i = 0; i < edges.Size(); i++ )
+	      sendarray. Append(edges[i] );
+
+	    for ( int i = 0; i < surfel.GetNP(); i++ )
+	      sendarray. Append(surfel[i] );
+	  }
+
+      }
+
+
+    if (id == 0)
+      MyMPI_Bcast ( sendarray );
+    else
+      MyMPI_Bcast ( recvarray );
+
+    cout << "have done bcast, id = " << id << endl;
+
+    if (id != 0)
+      {
+	Array<int,1> glob2loc_el;
+
 	glob2loc_el.SetSize (neglob);  
 	glob2loc_el = -1;
 	for ( int locel = 1; locel <= mesh.GetNE(); locel++)
 	  glob2loc_el[GetLoc2Glob_VolEl(locel)] = locel;
-      }
 
-
-    // MPI_Barrier (MPI_COMM_WORLD);
-
-    MPI_Request sendrequest;
-
-    if (id == 0)
-      {
- 	// PrintMessage (4, "UpdateCoarseGridGlobal : bcast, size = ", int (sendarray.Size()*sizeof(int)) );
-	MyMPI_Bcast ( sendarray );
-      }
-    else
-      sendrequest = MyMPI_ISend ( sendarray, 0, MPI_TAG_MESH  );
-
-
-    int nloops = (id == 0) ? ntasks-1 : 1;
-    for (int hi = 0; hi < nloops; hi++)
-      {
-	int sender;
-
-	if (id == 0)
-	  {
-	    sender = MyMPI_Recv ( recvarray, MPI_TAG_MESH );
-	    // PrintMessage (4, "have received from ", sender);
-	  }
-	else
-	  {
-	    MyMPI_Bcast ( recvarray );
-	    sender = 0;
-	  }
-	
-	// compare received vertices with own ones
 	int ii = 0;
-	// int cntel = 0;
-	int volel = 1;
-
-	if ( id != 0 )
-	  nfaglob = recvarray[ii++];
-	
-	if (id == 0) continue;
+	nfaglob = recvarray[ii++];
+	nedglob = recvarray[ii++];
 
 	Array<int> faces, edges;
 	Array<int> pnums, globalpnums;
@@ -632,16 +560,12 @@ namespace netgen
 	int recv_ne = recvarray[ii++];
 	for (int hi = 0; hi < recv_ne; hi++)
 	  {
-	    int globvolel   = recvarray[ii++];
+	    int globvolel = recvarray[ii++];
 	    int distnfa = recvarray[ii++];
 	    int distned = recvarray[ii++];
 	    int distnp  = recvarray[ii++];
 
-	    if ( id > 0 )      
-	      volel = glob2loc_el[globvolel];
-	    else
-	      volel = globvolel;
-	    
+	    int volel = glob2loc_el[globvolel];
 	    if (volel != -1)
 	      {
 		topology.GetElementFaces( volel, faces);
@@ -649,24 +573,67 @@ namespace netgen
 		const Element & volelement = mesh.VolumeElement (volel);
 
 		for ( int i = 0; i  < faces.Size(); i++)
-		  SetDistantFaceNum ( sender, faces[i], recvarray[ii++]);
+		  SetDistantFaceNum ( 0, faces[i], recvarray[ii++]);
 		
 		for ( int i = 0; i  < edges.Size(); i++)
-		  SetDistantEdgeNum ( sender, edges[i], recvarray[ii++]);
+		  SetDistantEdgeNum ( 0, edges[i], recvarray[ii++]);
 
-		for ( int i = 0; i  < distnp; i++)
-		  SetDistantPNum ( sender, volelement[i], recvarray[ii++]);
+		for ( int i = 0; i  < volelement.GetNP(); i++)
+		  SetDistantPNum ( 0, volelement[i], recvarray[ii++]);
 	      }
 	    else
 	      ii += distnfa + distned + distnp;
 	  }
+
+	
+	Array<int,1> glob2loc_sel;
+
+	int recv_nse = recvarray[ii++];
+	nseglob = recv_nse;
+
+	glob2loc_sel.SetSize (nseglob);  
+	glob2loc_sel = -1;
+	for ( int locel = 1; locel <= mesh.GetNSE(); locel++)
+	  glob2loc_sel[GetLoc2Glob_SurfEl(locel)] = locel;
+
+
+	for (int hi = 0; hi < recv_nse; hi++)
+	  {
+	    int globvolel = recvarray[ii++];
+	    int distned = recvarray[ii++];
+	    int distnp  = recvarray[ii++];
+
+	    int surfel = glob2loc_sel[globvolel];
+	    if (surfel != -1)
+	      {
+		topology.GetSurfaceElementEdges ( surfel, edges);
+		const Element2d & element = mesh.SurfaceElement (surfel);
+		
+		for ( int i = 0; i  < edges.Size(); i++)
+		  SetDistantEdgeNum ( 0, edges[i], recvarray[ii++]);
+
+		for ( int i = 0; i  < element.GetNP(); i++)
+		  SetDistantPNum ( 0, element[i], recvarray[ii++]);
+	      }
+	    else
+	      ii += distned + distnp;
+	  }
+
+
       }
-
-
-    coarseupdate = 1;
     
     if (id != 0)
-      MPI_Wait (&sendrequest, MPI_STATUS_IGNORE);
+      {
+	*testout << "l2d - vert = " << loc2distvert << endl;
+	*testout << "l2d - edge = " << loc2distedge << endl;
+	*testout << "l2d - el = " << loc2distel << endl;
+	*testout << "l2d - sel = " << loc2distsurfel << endl;
+      }
+
+    MPI_Barrier (MPI_COMM_WORLD);
+    cout << "update global complete, id = " << id << endl;
+    MPI_Barrier (MPI_COMM_WORLD);
+    coarseupdate = 1;
   }
 
 
@@ -676,6 +643,7 @@ namespace netgen
 
   void ParallelMeshTopology :: UpdateCoarseGrid ()
   {
+    cout << "updatecoarsegrid, id = " << id << endl;
     static int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid");
     NgProfiler::RegionTimer reg(timer);
 
@@ -684,6 +652,7 @@ namespace netgen
     if (id == 0)
       PrintMessage (1, "UPDATE COARSE GRID PARALLEL TOPOLOGY ");
 
+
     // find exchange edges - first send exchangeedges locnum, v1, v2
     // receive distant distnum, v1, v2
     // find matching
@@ -719,8 +688,10 @@ namespace netgen
     Array<int> cnt_send(ntasks-1);
 
 
-
-  
+    MPI_Barrier (MPI_HIGHORDER_COMM);
+    cout << "before ex edges" << endl;
+    *testout << "before ex edges" << endl;
+    MPI_Barrier (MPI_HIGHORDER_COMM);
 
     NgProfiler::StartTimer (timere);
 
@@ -793,6 +764,15 @@ namespace netgen
 
  
     NgProfiler::StopTimer (timere);
+
+    MPI_Barrier (MPI_HIGHORDER_COMM);
+    cout << "before ex faces" << endl;
+    *testout << "before ex faces" << endl;
+    MPI_Barrier (MPI_HIGHORDER_COMM);
+
+    if (mesh.GetDimension() == 3)
+      {
+
     NgProfiler::StartTimer (timerf);
 
     glob2loc.SetSize (nfaglob);
@@ -855,7 +835,12 @@ namespace netgen
 
 
     NgProfiler::StopTimer (timerf);
+      }
 
+    MPI_Barrier (MPI_HIGHORDER_COMM);
+    cout << "after ex faces" << endl;
+    *testout << "after ex faces" << endl;
+    MPI_Barrier (MPI_HIGHORDER_COMM);
 
 
     // set which elements are where for the master processor
diff --git a/libsrc/meshing/paralleltop.hpp b/libsrc/meshing/paralleltop.hpp
index 65dd67c6..52c03c52 100644
--- a/libsrc/meshing/paralleltop.hpp
+++ b/libsrc/meshing/paralleltop.hpp
@@ -17,9 +17,10 @@ namespace netgen
     int nseg, nsurfel;
 
     // number of global elements, vertices, ???,  faces
-    int neglob, nvglob;
+    int neglob, nseglob, nvglob;
     int nparel;
     int nfaglob;
+    int nedglob;
 
     /**
        mapping from local to distant vertex number
@@ -45,6 +46,7 @@ namespace netgen
 
     void SetNVGlob ( int anvglob )   { nvglob = anvglob; }
     void SetNEGlob ( int aneglob )   { neglob = aneglob; }
+    void SetNSEGlob ( int anseglob )   { nseglob = anseglob; }
 
     int GetNVGlob ()  { return nvglob; }
     int GetNEGlob ()  { return neglob; }
@@ -59,6 +61,7 @@ namespace netgen
 
     int GetLoc2Glob_Vert  ( int locnum ) const { return loc2distvert[locnum][0]; }
     int GetLoc2Glob_VolEl ( int locnum ) const { return loc2distel[locnum-1][0]; }
+    int GetLoc2Glob_SurfEl ( int locnum ) const { return loc2distsurfel[locnum-1][0]; }
 
     void GetVertNeighbours ( int vnum, Array<int> & dests ) const;
 
diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp
index 5f46d553..7e638bcd 100644
--- a/libsrc/meshing/topology.cpp
+++ b/libsrc/meshing/topology.cpp
@@ -1114,50 +1114,10 @@ namespace netgen
 		    if ( ntasks > 1 )
 		      {
 			if ( !paralleltop.DoCoarseUpdate() ) continue;
-			// "illegal" faces are exchange faces
-			/*
-			  (*testout) << "exchange face : " << i << endl;
-			  (*testout) << "points = " << face2vert[i] << endl;
-			*/
-			// 		  (*testout) << "global points = ";
-			// 		  for ( int j = 0; j < 3; j++ )
-			// // 		    (*testout) << face2vert[i].I(j+1) << " -> " 
-			// // 			       << paralleltop.GetLoc2Glob_Vert( face2vert[i].I(j+1) ) << ",  ";
-			// 		  (*testout) << endl;
-			// 		  if ( !paralleltop.IsExchangeFace (i+1) )
-			// 		    paralleltop.SetRefinementFace (i+1);
-
-
-			/*
-			paralleltop.SetExchangeFace (i+1);
-		  
-			for (int j = 0; j < 4; j++)		    
-			  {
-			    if ( face2vert[i].I(j+1) > 0 )
-			      paralleltop.SetExchangeVert(face2vert[i].I(j+1));
-			  }
-		  
-			Array<int> faceedges;
-			GetFaceEdges (i+1, faceedges);
-			for ( int j = 0; j < faceedges.Size(); j++)
-			  {
-			    paralleltop.SetExchangeEdge ( faceedges[j] );
-			    int v1, v2;
-			    GetEdgeVertices(faceedges[j], v1, v2 );
-			  }
-			*/
-
-			/*
-			  (*testout) << "pos = ";
-			  for (int j = 0; j < 4; j++)
-			  if (face2vert[i].I(j+1) >= 1)
-			  (*testout) << mesh[(PointIndex)face2vert[i].I(j+1)] << " ";
-			  (*testout) << endl;
-			*/
 		      }
 		    else
-		      {
 #endif
+		      {
 			(*testout) << "illegal face : " << i << endl;
 			(*testout) << "points = " << face2vert[i] << endl;
 			(*testout) << "pos = ";
@@ -1185,42 +1145,36 @@ namespace netgen
 
 				}
 			  }
-#ifdef PARALLEL
 		      }
-#endif
 		  }
 	      }
 
-#ifndef PARALLEL
-	    if (cnt_err)
-	      cout << cnt_err << " elements are not matching !!!" << endl;
-#else
 	    if (cnt_err && ntasks == 1)
 	      cout << cnt_err << " elements are not matching !!!" << endl;
-	    else if (cnt_err && ntasks > 1)
-	      {
-		// cout << "p" << id << ":  " << cnt_err << " elements are not local" << endl;
-		isparallel = 1;
-	      }
-	    else if ( ntasks > 1 )
-	      ; 
-	    // cout << "p" << id << ":  " << "Partition " << id << " is totally local" << endl;
-#endif
+
+	    if (cnt_err && ntasks > 1)
+	      isparallel = 1;
 	  }
+      }
+    
 
 #ifdef PARALLEL
-	if ( isparallel )
-	  {
-	    paralleltop.Update();
-	    if ( paralleltop.DoCoarseUpdate() )
-	      paralleltop.UpdateCoarseGrid();
-	  }
-#endif
-
-
-
-
+    cout << "is par = " << int(isparallel) << ", id = " << id << endl;
+    if (mesh.GetDimension() == 3)
+      if (isparallel != (id != 0))
+	{
+	  cerr << " ****************************** " << endl;
+	  cerr << " **** wrong isparallel   ****** " << endl;
+	  cerr << " ****************************** " << endl;
+	}
+    
+    if (id != 0)    // if (isparallel)
+      {
+	paralleltop.Update();
+	if ( paralleltop.DoCoarseUpdate() )
+	  paralleltop.UpdateCoarseGrid();
       }
+#endif
  
  
   
@@ -1737,7 +1691,8 @@ namespace netgen
 
 
 
-  void MeshTopology :: GetSegmentVolumeElements ( int segnr, Array<int> & volels ) const
+  void MeshTopology :: 
+  GetSegmentVolumeElements ( int segnr, Array<int> & volels ) const
   {
     int v1, v2;
     GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 );
@@ -1749,6 +1704,24 @@ namespace netgen
     for ( int eli1=1; eli1 <= volels1.Size(); eli1++)
       if ( volels2.Contains( volels1.Elem(eli1) ) )
 	volels.Append ( volels1.Elem(eli1) );
-
   }
+
+  void MeshTopology :: 
+  GetSegmentSurfaceElements (int segnr, Array<int> & els) const
+  {
+    int v1, v2;
+    GetEdgeVertices ( GetSegmentEdge (segnr), v1, v2 );
+    Array<int> els1, els2;
+    GetVertexSurfaceElements ( v1, els1 );
+    GetVertexSurfaceElements ( v2, els2 );
+    els.SetSize(0);
+
+    for ( int eli1=1; eli1 <= els1.Size(); eli1++)
+      if ( els2.Contains( els1.Elem(eli1) ) )
+	els.Append ( els1.Elem(eli1) );
+  }
+
+
+
+
 }
diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp
index 92e7c182..a453fa84 100644
--- a/libsrc/meshing/topology.hpp
+++ b/libsrc/meshing/topology.hpp
@@ -124,7 +124,8 @@ public:
 
   
   int GetVerticesEdge ( int v1, int v2) const;
-  void GetSegmentVolumeElements ( int segnr, Array<int> & surfels ) const;
+  void GetSegmentVolumeElements ( int segnr, Array<int> & els ) const;
+  void GetSegmentSurfaceElements ( int segnr, Array<int> & els ) const;
 };
 
 
diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp
index a8eba0df..5c5f61f9 100644
--- a/libsrc/visualization/vssolution.cpp
+++ b/libsrc/visualization/vssolution.cpp
@@ -2192,7 +2192,7 @@ namespace netgen
     if (minv == maxv) maxv = minv+1e-6;
 
 #ifdef PARALLEL
-    if (id == 0)
+    if ((ntasks > 1) && (id == 0))
       {
 	minv = 1e99;
 	maxv = -1e99;