diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp
index 8947e879..eeb2a083 100644
--- a/libsrc/meshing/bisect.cpp
+++ b/libsrc/meshing/bisect.cpp
@@ -3708,6 +3708,7 @@ namespace netgen
     mesh.SurfaceElements().SetAllocSize (mtris.Size()+mquads.Size());    
     NgProfiler::StopTimer (timer2a);        
     NgProfiler::StartTimer (timer2b);
+    /*
     for (auto & trig : mtris)
       {
 	Element2d el(TRIG);
@@ -3720,6 +3721,28 @@ namespace netgen
 	  }
 	mesh.AddSurfaceElement (el);
       }
+    */
+
+    mesh.SurfaceElements().SetSize(mtris.Size());
+    // for (size_t i = 0; i < mtris.Size(); i++)
+    ParallelForRange
+      (opt.task_manager, mtris.Size(), [&] (size_t begin, size_t end)
+       {
+         for (size_t i = begin; i < end; i++)
+          {
+            Element2d el(TRIG);
+            auto & trig = mtris[i];
+            el.SetIndex (trig.surfid);
+            el.SetOrder (trig.order);
+            for (int j = 0; j < 3; j++)
+              {
+                el[j] = trig.pnums[j];
+                el.GeomInfoPi(j+1) = trig.pgeominfo[j];
+              }
+            mesh.SetSurfaceElement (SurfaceElementIndex(i), el);
+          }
+       });
+    
     for (int i = 1; i <= mquads.Size(); i++)
       {
 	Element2d el(QUAD);
diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp
index d444a8b1..1db47a41 100644
--- a/libsrc/meshing/meshclass.cpp
+++ b/libsrc/meshing/meshclass.cpp
@@ -325,6 +325,42 @@ namespace netgen
     return si;
   }
 
+  void Mesh :: SetSurfaceElement (SurfaceElementIndex sei, const Element2d & el)
+  {
+    int maxn = el[0];
+    for (int i = 1; i < el.GetNP(); i++)
+      if (el[i] > maxn) maxn = el[i];
+
+    maxn += 1-PointIndex::BASE;
+
+    if (maxn <= points.Size())
+      {
+        for (int i = 0; i < el.GetNP(); i++)
+          if (points[el[i]].Type() > SURFACEPOINT)
+            points[el[i]].SetType(SURFACEPOINT);
+      }
+
+    surfelements[sei] = el;
+    if (el.index > facedecoding.Size())
+      cerr << "has no facedecoding: fd.size = " << facedecoding.Size() << ", ind = " << el.index << endl;
+
+    /*
+    surfelements.Last().next = facedecoding[el.index-1].firstelement;
+    facedecoding[el.index-1].firstelement = sei;
+    */
+
+    // add lock-free to list
+    surfelements[sei].next = facedecoding[el.index-1].firstelement;
+    auto & head = reinterpret_cast<atomic<SurfaceElementIndex>&> (facedecoding[el.index-1].firstelement);
+    while (!head.compare_exchange_weak (surfelements[sei].next, sei))
+      ;
+    
+    /*
+    if (SurfaceArea().Valid())
+      SurfaceArea().Add (el);
+    */
+  }
+
 
   ElementIndex Mesh :: AddVolumeElement (const Element & el)
   { 
diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp
index 89834d71..162dfd79 100644
--- a/libsrc/meshing/meshclass.hpp
+++ b/libsrc/meshing/meshclass.hpp
@@ -261,6 +261,8 @@ namespace netgen
     Array<Element0d> pointelements;  // only via python interface
 
     DLL_HEADER SurfaceElementIndex AddSurfaceElement (const Element2d & el);
+    // write to pre-allocated container, thread-safe
+    DLL_HEADER void SetSurfaceElement (SurfaceElementIndex sei, const Element2d & el);
     
     // [[deprecated("Use DeleteSurfaceElement(SurfaceElementIndex) instead of int !")]]
     void DeleteSurfaceElement (int eli)
diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp
index 3fe659bc..2eb2e65a 100644
--- a/libsrc/meshing/meshtype.hpp
+++ b/libsrc/meshing/meshtype.hpp
@@ -228,8 +228,11 @@ namespace netgen
   public:
     SurfaceElementIndex () { ; }
     SurfaceElementIndex (int ai) : i(ai) { ; }
+    /*
     SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) 
     { i = ai.i; return *this; }
+    */
+    SurfaceElementIndex & operator= (const SurfaceElementIndex & ai) = default;
     SurfaceElementIndex & operator= (int ai) { i = ai; return *this; }
     operator int () const { return i; }
     SurfaceElementIndex operator++ (int) { SurfaceElementIndex hi(*this); i++; return hi; }