diff --git a/libsrc/meshing/clusters.cpp b/libsrc/meshing/clusters.cpp
index 0c2dd7bb..52e5eec1 100644
--- a/libsrc/meshing/clusters.cpp
+++ b/libsrc/meshing/clusters.cpp
@@ -220,24 +220,7 @@ namespace netgen
 	  {
 	    const Element & el = mesh.VolumeElement(i);
 	    ELEMENT_TYPE typ = el.GetType();
-	  
-	    top.GetElementEdges (i, ednums);
-	    top.GetElementFaces (i, fanums);
-	  
-	    int elnv = top.GetNVertices (typ);
-	    int elned = ednums.Size();
-	    int elnfa = fanums.Size();
-	  
-	    nnums.SetSize(elnv+elned+elnfa+1);
-	    for (int j = 1; j <= elnv; j++)
-	      nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
-	    for (int j = 1; j <= elned; j++)
-	      nnums.Elem(elnv+j) = nv+ednums.Elem(j);
-	    for (int j = 1; j <= elnfa; j++)
-	      nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j);
-	    nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i;
-
-	  
+	  	  
 	    const int * clustertab = NULL;
 	    switch (typ)
 	      {
@@ -280,6 +263,25 @@ namespace netgen
 	      }
 	  
 	    if (clustertab)
+              {
+                top.GetElementEdges (i, ednums);
+                top.GetElementFaces (i, fanums);
+                
+                int elnv = top.GetNVertices (typ);
+                int elned = ednums.Size();
+                int elnfa = fanums.Size();
+                
+                nnums.SetSize(elnv+elned+elnfa+1);
+                for (int j = 1; j <= elnv; j++)
+                  nnums.Elem(j) = el.PNum(j)+1-PointIndex::BASE;
+                for (int j = 1; j <= elned; j++)
+                  nnums.Elem(elnv+j) = nv+ednums.Elem(j);
+                for (int j = 1; j <= elnfa; j++)
+                  nnums.Elem(elnv+elned+j) = nv+ned+fanums.Elem(j);
+                nnums.Elem(elnv+elned+elnfa+1) = nv+ned+nfa+i;
+                
+
+                
 	      for (int j = 0; j < nnums.Size(); j++)
 		for (int k = 0; k < j; k++)
 		  if (clustertab[j] == clustertab[k])
@@ -312,6 +314,7 @@ namespace netgen
 			  changed = 1;
 			}
 		    }
+              }
 
 	    /*
 	      if (clustertab)
diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp
index 752d16a0..2fd6b7c9 100644
--- a/libsrc/meshing/topology.cpp
+++ b/libsrc/meshing/topology.cpp
@@ -1096,6 +1096,7 @@ namespace netgen
 	    surf2volelement.Elem(i)[0] = 0;
 	    surf2volelement.Elem(i)[1] = 0;
 	  }
+        (*tracer) ("Topology::Update build surf2vol", false);        
 	for (int i = 1; i <= ne; i++)
 	  for (int j = 0; j < 6; j++)
 	    {
@@ -1109,6 +1110,7 @@ namespace netgen
 		  surf2volelement.Elem(sel)[0] = i;
 		}
 	    }
+        (*tracer) ("Topology::Update build surf2vol", true);        
 
 	face2vert.SetAllocSize (face2vert.Size());
 
@@ -1120,9 +1122,11 @@ namespace netgen
 	// paralleltop.Reset ();
 #endif
 
+        (*tracer) ("Topology::Update count face_els", false);
 	Array<short int> face_els(nfa), face_surfels(nfa);
 	face_els = 0;
 	face_surfels = 0;
+        /*
 	Array<int> hfaces;
 	for (int i = 1; i <= ne; i++)
 	  {
@@ -1130,8 +1134,22 @@ namespace netgen
 	    for (int j = 0; j < hfaces.Size(); j++)
 	      face_els[hfaces[j]-1]++;
 	  }
+        */
+        ParallelForRange
+          (tm, ne,
+            [&] (size_t begin, size_t end)
+            {
+              Array<int> hfaces;              
+              for (ElementIndex ei = begin; ei < end; ei++)
+                {
+                  GetElementFaces (ei+1, hfaces);
+                  for (auto f : hfaces)
+                    AsAtomic(face_els[f-1])++;
+                }
+            });        
 	for (int i = 1; i <= nse; i++)
 	  face_surfels[GetSurfaceElementFace (i)-1]++;
+        (*tracer) ("Topology::Update count face_els", true);
 
 
 	if (ne)