Merge branch 'parallel_meshing' into 'master'

Parallel meshing

See merge request jschoeberl/netgen!398

libsrc/core/taskmanager.cpp

@@ -235,13 +235,14 @@ namespace ngcore
     const function<void(TaskInfo&)> * func;
     int mynr;
     int total;
+    int producing_thread;
     atomic<int> * endcnt;
 
     TNestedTask () { ; }
     TNestedTask (const function<void(TaskInfo&)> & _func,
                  int _mynr, int _total,
-                 atomic<int> & _endcnt)
+                 atomic<int> & _endcnt, int prod_tid)
-      : func(&_func), mynr(_mynr), total(_total), endcnt(&_endcnt)
+      : func(&_func), mynr(_mynr), total(_total), endcnt(&_endcnt), producing_thread(prod_tid)
     {
       ;
     }
@@ -260,12 +261,14 @@ namespace ngcore
     TPToken ptoken(taskqueue); 
 
     int num = endcnt;
+    auto tid = TaskManager::GetThreadId();
     for (int i = 0; i < num; i++)
-      taskqueue.enqueue (ptoken, { afunc, i, num, endcnt });
+      taskqueue.enqueue (ptoken, { afunc, i, num, endcnt, tid });
   }
 
   bool TaskManager :: ProcessTask()
   {
+    static Timer t("process task");
     TNestedTask task;
     TCToken ctoken(taskqueue); 
     
@@ -282,8 +285,14 @@ namespace ngcore
           cout << "process nested, nr = " << ti.task_nr << "/" << ti.ntasks << endl;
         }
         */
+        if(trace && task.producing_thread != ti.thread_nr)
+            trace->StartTask (ti.thread_nr, t, PajeTrace::Task::ID_TIMER, task.producing_thread);
+
         (*task.func)(ti);
         --*task.endcnt;
+
+        if(trace && task.producing_thread != ti.thread_nr)
+            trace->StopTask (ti.thread_nr, t);
         return true;
       }
     return false;

libsrc/meshing/debugging.hpp

@@ -0,0 +1,51 @@
+#include <meshclass.hpp>
+
+
+namespace netgen
+{
+    inline unique_ptr<Mesh> GetOpenElements( const Mesh & m, int dom = 0 )
+    {
+        static Timer t("GetOpenElements"); RegionTimer rt(t);
+        auto mesh = make_unique<Mesh>();
+        *mesh = m;
+
+        Array<bool, PointIndex> interesting_points(mesh->GetNP());
+        interesting_points = false;
+
+        mesh->FindOpenElements(dom);
+        NgArray<Element2d> openelements;
+        openelements = mesh->OpenElements();
+
+        for (auto & el : openelements)
+            for (auto i : el.PNums())
+                interesting_points[i] = true;
+
+        for (auto & el : mesh->VolumeElements())
+        {
+            int num_interesting_points = 0;
+
+            for (auto pi : el.PNums())
+                if(interesting_points[pi])
+                    num_interesting_points++;
+
+            if(num_interesting_points==0)
+                el.Delete();
+            el.SetIndex(num_interesting_points);
+        }
+
+        mesh->SetMaterial(1, "1_point");
+        mesh->SetMaterial(2, "2_points");
+        mesh->SetMaterial(3, "3_points");
+        mesh->SetMaterial(4, "4_points");
+        mesh->Compress();
+
+        mesh->ClearSurfaceElements();
+
+        for (auto & el : openelements)
+            mesh->AddSurfaceElement( el );
+
+        return mesh;
+    }
+
+
+}

libsrc/meshing/delaunay.cpp

@@ -235,9 +235,11 @@ namespace netgen
 			 NgArray<int> & freelist, SphereList & list,
 			 IndexSet & insphere, IndexSet & closesphere, Array<DelaunayTet> & newels)
   {
-    static Timer t("Meshing3::AddDelaunayPoint");// RegionTimer reg(t);
+    static Timer t("Meshing3::AddDelaunayPoint", NoTracing, NoTiming); RegionTimer reg(t);
-    // static Timer tsearch("addpoint, search");
+    static Timer tsearch("addpoint, search", NoTracing, NoTiming);
-    // static Timer tinsert("addpoint, insert");
+    static Timer tfind("addpoint, find all tets", NoTracing, NoTiming);
+    static Timer tnewtets("addpoint, build new tets", NoTracing, NoTiming);
+    static Timer tinsert("addpoint, insert", NoTracing, NoTiming);
 
     /*
       find any sphere, such that newp is contained in
@@ -277,7 +279,7 @@ namespace netgen
       }
     */
 
-    // tsearch.Start();
+    tsearch.Start();
     double minquot{1e20};
     tettree.GetFirstIntersecting
       (newp, newp, [&](const auto pi)
@@ -300,7 +302,7 @@ namespace netgen
           }
         return false;
        } );
-    // tsearch.Stop();
+    tsearch.Stop();
 
     if (cfelind == -1)
       {
@@ -308,6 +310,7 @@ namespace netgen
 	return;
       }
 	
+    tfind.Start();
     /*
       insphere:     point is in sphere -> delete element
       closesphere:  point is close to sphere -> considered for same center
@@ -399,6 +402,8 @@ namespace netgen
 	    }
       } // while (changed)
 
+    tfind.Stop();
+    tnewtets.Start();
     newels.SetSize(0);
     
     Element2d face(TRIG);
@@ -553,10 +558,11 @@ namespace netgen
 	    tpmax.SetToMax (*pp[k]);
 	  }
 	tpmax = tpmax + 0.01 * (tpmax - tpmin);
-        // tinsert.Start();
+        tinsert.Start();
 	tettree.Insert (tpmin, tpmax, nelind);
-        // tinsert.Stop();
+        tinsert.Stop();
       }
+      tnewtets.Stop();
   }
 
 
@@ -1539,7 +1545,7 @@ namespace netgen
 
     PrintMessage (1, "Delaunay meshing");
     PrintMessage (3, "number of points: ", mesh.GetNP());
-    PushStatus ("Delaunay meshing");
+    // PushStatus ("Delaunay meshing");
 
 
     NgArray<DelaunayTet> tempels;
@@ -1627,20 +1633,20 @@ namespace netgen
       // tempmesh.Save ("tempmesh.vol");
 
       {
-        RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
+        MeshOptimize3d meshopt(mp);
-        for (int i = 1; i <= 4; i++)
+        tempmesh.Compress();
-          {
         tempmesh.FindOpenElements ();
-
+        RegionTaskManager rtm(mp.parallel_meshing ? mp.nthreads : 0);
+        for (auto i : Range(10))
+          {
             PrintMessage (5, "Num open: ", tempmesh.GetNOpenElements());
-            tempmesh.CalcSurfacesOfNode ();
 
+            if(i%5==0)
                 tempmesh.FreeOpenElementsEnvironment (1);
 
-            MeshOptimize3d meshopt(mp);
-            // tempmesh.CalcSurfacesOfNode();
             meshopt.SwapImprove(tempmesh, OPT_CONFORM);
           }
+        tempmesh.Compress();
       }
     
       MeshQuality3d (tempmesh);
@@ -1669,6 +1675,6 @@ namespace netgen
     mesh.FindOpenElements(domainnr);
 
     mesh.Compress();
-    PopStatus ();
+    // PopStatus ();
   }
 }

libsrc/meshing/improve2.hpp

@@ -1,15 +1,42 @@
 #ifndef FILE_IMPROVE2
 #define FILE_IMPROVE2
 
+inline void AppendEdges( const Element2d & elem, PointIndex pi, Array<std::tuple<PointIndex,PointIndex>> & edges )
+{
+  for (int j = 0; j < 3; j++)
+  {
+      PointIndex pi0 = elem[j];
+      PointIndex pi1 = elem[(j+1)%3];
+      if (pi1 < pi0) Swap(pi0, pi1);
+      if(pi0==pi)
+          edges.Append(std::make_tuple(pi0, pi1));
+  }
+}
+
+inline void AppendEdges( const Element & elem, PointIndex pi, Array<std::tuple<PointIndex,PointIndex>> & edges )
+{
+  static constexpr int tetedges[6][2] =
+  { { 0, 1 }, { 0, 2 }, { 0, 3 },
+      { 1, 2 }, { 1, 3 }, { 2, 3 } };
+
+  if(elem.flags.fixed)
+      return;
+  for (int j = 0; j < 6; j++)
+  {
+      PointIndex pi0 = elem[tetedges[j][0]];
+      PointIndex pi1 = elem[tetedges[j][1]];
+      if (pi1 < pi0) Swap(pi0, pi1);
+      if(pi0==pi)
+          edges.Append(std::make_tuple(pi0, pi1));
+  }
+}
+
 template<typename TINDEX>
 void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & elementsonnode, Array<std::tuple<PointIndex, PointIndex>> & edges )
 {
+  static_assert(is_same_v<TINDEX, ElementIndex>||is_same_v<TINDEX,SurfaceElementIndex>, "Invalid type for TINDEX");
   static Timer tbuild_edges("Build edges"); RegionTimer reg(tbuild_edges);
 
-  static constexpr int tetedges[6][2] =
-    { { 0, 1 }, { 0, 2 }, { 0, 3 },
-        { 1, 2 }, { 1, 3 }, { 2, 3 } };
-
   int ntasks = 4*ngcore::TaskManager::GetMaxThreads();
   Array<Array<std::tuple<PointIndex,PointIndex>>> task_edges(ntasks);
 
@@ -26,29 +53,7 @@ void BuildEdgeList( const Mesh & mesh, const Table<TINDEX, PointIndex> & element
             const auto & elem = mesh[ei];
             if (elem.IsDeleted()) continue;
 
-            static_assert(is_same_v<TINDEX, ElementIndex>||is_same_v<TINDEX,SurfaceElementIndex>, "Invalid type for TINDEX");
+            AppendEdges(elem, pi, local_edges);
-            if constexpr(is_same_v<TINDEX, SurfaceElementIndex>)
-            {
-                for (int j = 0; j < 3; j++)
-                {
-                    PointIndex pi0 = elem[j];
-                    PointIndex pi1 = elem[(j+1)%3];
-                    if (pi1 < pi0) Swap(pi0, pi1);
-                    if(pi0==pi)
-                        local_edges.Append(std::make_tuple(pi0, pi1));
-                }
-            }
-            else if constexpr(is_same_v<TINDEX, ElementIndex>)
-            {
-                for (int j = 0; j < 6; j++)
-                {
-                    PointIndex pi0 = elem[tetedges[j][0]];
-                    PointIndex pi1 = elem[tetedges[j][1]];
-                    if (pi1 < pi0) Swap(pi0, pi1);
-                    if(pi0==pi)
-                        local_edges.Append(std::make_tuple(pi0, pi1));
-                }
-            }
         }
         QuickSort(local_edges);
 

libsrc/meshing/improve3.cpp

@@ -2717,8 +2717,24 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
   int ne = mesh.GetNE();
 
   mesh.BuildBoundaryEdges(false);
+  BitArray free_points(mesh.GetNP()+PointIndex::BASE);
+  free_points.Clear();
 
-  auto elementsonnode = mesh.CreatePoint2ElementTable();
+  ParallelForRange(mesh.VolumeElements().Range(), [&] (auto myrange)
+      {
+        for (ElementIndex eli : myrange)
+          {
+            const auto & el = mesh[eli];
+            if(el.flags.fixed)
+              continue;
+
+            for (auto pi : el.PNums())
+              if(!free_points[pi])
+                  free_points.SetBitAtomic(pi);
+          }
+      });
+
+  auto elementsonnode = mesh.CreatePoint2ElementTable(free_points);
 
   NgArray<ElementIndex> hasbothpoints;
 
@@ -2736,7 +2752,7 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
 	  const Element2d & hel = mesh.OpenElement(i);
 	  INDEX_3 face(hel[0], hel[1], hel[2]);
 	  face.Sort();
-	  faces.Set (face, 1);
+	  faces.Set (face, i);
 	}
     }
 
@@ -2755,6 +2771,8 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
 
   tloop.Start();
 
+  auto num_elements_before = mesh.VolumeElements().Range().Next();
+
   ParallelForRange(Range(edges), [&] (auto myrange)
   {
     for(auto i : myrange)
@@ -2786,6 +2804,30 @@ void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal,
 
   PrintMessage (5, cnt, " swaps performed");
 
+  if(goal == OPT_CONFORM)
+  {
+      // Remove open elements that were closed by new tets
+      auto & open_els = mesh.OpenElements();
+
+      for (auto & el : mesh.VolumeElements().Range( num_elements_before, mesh.VolumeElements().Range().Next() ))
+      {
+          for (auto i : Range(1,5))
+          {
+              Element2d sel;
+              el.GetFace(i, sel);
+              INDEX_3 face(sel[0], sel[1], sel[2]);
+              face.Sort();
+              if(faces.Used(face))
+                  open_els[faces.Get(face)-1].Delete();
+          }
+      }
+
+      for(int i=open_els.Size()-1; i>=0; i--)
+          if(open_els[i].IsDeleted())
+              open_els.Delete(i);
+
+      mesh.DeleteBoundaryEdges();
+  }
   mesh.Compress ();
 
   multithread.task = savetask;

libsrc/meshing/localh.cpp

@@ -91,6 +91,93 @@ namespace netgen
     delete root;
   }
 
+  unique_ptr<LocalH> LocalH :: Copy ()
+  {
+    static Timer t("LocalH::Copy"); RegionTimer rt(t);
+    auto lh = make_unique<LocalH>(boundingbox, grading, dimension);
+    std::map<GradingBox*, GradingBox*> mapping;
+    lh->boxes.SetSize(boxes.Size());
+
+    for(auto i : boxes.Range())
+    {
+      lh->boxes[i] = new GradingBox();
+      auto & bnew = *lh->boxes[i];
+      auto & b = *boxes[i];
+      bnew.xmid[0] = b.xmid[0];
+      bnew.xmid[1] = b.xmid[1];
+      bnew.xmid[2] = b.xmid[2];
+      bnew.h2 = b.h2;
+      bnew.hopt = b.hopt;
+      bnew.flags = b.flags;
+      mapping[&b] = &bnew;
+    }
+
+    for(auto i : boxes.Range())
+    {
+      auto & bnew = *lh->boxes[i];
+      auto & b = *boxes[i];
+      for(auto k : Range(8))
+      {
+        if(b.childs[k])
+          bnew.childs[k] = mapping[b.childs[k]];
+      }
+
+      if(b.father)
+        bnew.father = mapping[b.father];
+    }
+
+    lh->root = mapping[root];
+    return lh;
+  }
+
+  unique_ptr<LocalH> LocalH :: Copy( const Box<3> & bbox )
+  {
+    static Timer t("LocalH::Copy with bounding box"); RegionTimer rt(t);
+    auto lh = make_unique<LocalH>(boundingbox, grading, dimension);
+    std::map<GradingBox*, GradingBox*> mapping;
+    lh->boxes.SetAllocSize(boxes.Size());
+
+    for(auto i : boxes.Range())
+    {
+      auto & b = *boxes[i];
+      auto h = b.H2();
+      Vec<3> vh = {h,h,h};
+      Box<3> box( b.PMid() - vh, b.PMid() + vh);
+      if(!box.Intersect(bbox))
+          continue;
+      lh->boxes.Append(new GradingBox());
+      auto & bnew = *lh->boxes.Last();
+      bnew.xmid[0] = b.xmid[0];
+      bnew.xmid[1] = b.xmid[1];
+      bnew.xmid[2] = b.xmid[2];
+      bnew.h2 = b.h2;
+      bnew.hopt = b.hopt;
+      bnew.flags = b.flags;
+      mapping[&b] = &bnew;
+    }
+
+    for(auto i : boxes.Range())
+    {
+      auto & b = *boxes[i];
+      if(mapping.count(&b)==0)
+          continue;
+
+      auto & bnew = *mapping[&b];
+      for(auto k : Range(8))
+      {
+        if(b.childs[k] && mapping.count(b.childs[k]))
+          bnew.childs[k] = mapping[b.childs[k]];
+      }
+
+      if(b.father && mapping.count(b.father))
+        bnew.father = mapping[b.father];
+    }
+
+    lh->root = mapping[root];
+    return lh;
+
+  }
+
   void LocalH :: Delete ()
   {
     root->DeleteChilds();
@@ -405,8 +492,8 @@ namespace netgen
   void LocalH :: FindInnerBoxes (AdFront3 * adfront,
 				 int (*testinner)(const Point3d & p1))
   {
-    static int timer = NgProfiler::CreateTimer ("LocalH::FindInnerBoxes");
+    static Timer timer("LocalH::FindInnerBoxes");
-    NgProfiler::RegionTimer reg (timer);
+    RegionTimer reg (timer);
 
 
     int nf = adfront->GetNF();
@@ -812,6 +899,7 @@ namespace netgen
   
   void LocalH :: GetOuterPoints (NgArray<Point<3> > & points)
   {
+    static Timer t("LocalH::GetOuterPoints"); RegionTimer rt(t);
     for (int i = 0; i < boxes.Size(); i++)
       if (!boxes[i]->flags.isinner && !boxes[i]->flags.cutboundary)
 	points.Append ( boxes[i] -> PMid());

libsrc/meshing/localh.hpp

@@ -98,6 +98,9 @@ namespace netgen
 
     ~LocalH();
     ///
+    unique_ptr<LocalH> Copy();
+    unique_ptr<LocalH> Copy( const Box<3> & bbox );
+    ///
     void Delete();
     ///
     void DoArchive(Archive& ar);
@@ -189,6 +192,8 @@ namespace netgen
     ///
     void ConvexifyRec (GradingBox * box);
 
+    unique_ptr<LocalH> CopyRec( const Box<3> & bbox, GradingBox * current );
+
     friend ostream & operator<< (ostream & ost, const LocalH & loch);
   };
 

libsrc/meshing/meshclass.cpp

@@ -384,6 +384,8 @@ namespace netgen
         volelements.Append (el);
       }
     volelements.Last().flags.illegal_valid = 0;
+    volelements.Last().flags.fixed = 0;
+    volelements.Last().flags.deleted = 0;
 
     // while (volelements.Size() > eltyps.Size())
     // eltyps.Append (FREEELEMENT);
@@ -405,6 +407,8 @@ namespace netgen
 
     volelements[ei]  = el;
     volelements[ei].flags.illegal_valid = 0;
+    volelements[ei].flags.fixed = 0;
+    volelements[ei].flags.deleted = 0;
   }
 
 
@@ -436,6 +440,7 @@ namespace netgen
 
   void Mesh :: Save (ostream & outfile) const
   {
+    static Timer timer("Mesh::Save"); RegionTimer rt(timer);
     int i, j;
 
     double scale = 1;  // globflags.GetNumFlag ("scale", 1);
@@ -2910,6 +2915,7 @@ namespace netgen
 
   void Mesh :: FreeOpenElementsEnvironment (int layers)
   {
+    static Timer timer("FreeOpenElementsEnvironment"); RegionTimer rt(timer);
     int i, j, k;
     PointIndex pi;
     const int large = 9999;
@@ -6570,6 +6576,8 @@ namespace netgen
                [&](auto & table, ElementIndex ei)
                {
                  const auto & el = (*this)[ei];
+                 if(el.IsDeleted())
+                     return;
 
                  for (PointIndex pi : el.PNums())
                    if(free_points[pi])
@@ -6581,6 +6589,8 @@ namespace netgen
                [&](auto & table, ElementIndex ei)
                {
                  const auto & el = (*this)[ei];
+                 if(el.IsDeleted())
+                     return;
 
                  for (PointIndex pi : el.PNums())
                    table.Add (pi, ei);

libsrc/meshing/meshclass.hpp

@@ -487,6 +487,8 @@ namespace netgen
 
     auto & OpenElements() const { return openelements; }
 
+    auto & OpenElements() { return openelements; }
+    
     /// are also quads open elements
     bool HasOpenQuads () const;
 
@@ -510,6 +512,11 @@ namespace netgen
       return boundaryedges->Used (i2);
     }
 
+    void DeleteBoundaryEdges ()
+    {
+        boundaryedges = nullptr;
+    }
+
     bool IsSegment (PointIndex pi1, PointIndex pi2) const
     {
       INDEX_2 i2 (pi1, pi2);

libsrc/meshing/meshfunc.cpp

@@ -1,5 +1,6 @@
 #include <mystdlib.h>
 #include "meshing.hpp"
+#include "debugging.hpp"
 
 namespace netgen
 {
@@ -10,106 +11,182 @@ namespace netgen
   extern const char * pyramidrules2[];
   extern const char * hexrules[];
 
-
+  struct MeshingData
-  // extern double teterrpow; 
-  MESHING3_RESULT MeshVolume (const MeshingParameters & c_mp, Mesh& mesh3d)
   {
-    static Timer t("MeshVolume"); RegionTimer reg(t);
+      int domain;
+
+      // mesh for one domain (contains all adjacent surface elments)
+      unique_ptr<Mesh> mesh;
+
+      // maps from local (domain) mesh to global mesh
+      Array<PointIndex, PointIndex> pmap;
+
+      // Array<INDEX_2> connected_pairs;
+
+      MeshingParameters mp;
+
+      unique_ptr<Meshing3> meshing;
+  };
+
+  // extract surface meshes belonging to individual domains
+  Array<MeshingData> DivideMesh(Mesh & mesh, const MeshingParameters & mp)
+  {
+      static Timer timer("DivideMesh"); RegionTimer rt(timer);
+
+      Array<MeshingData> ret;
+      auto num_domains = mesh.GetNDomains();
+
+      if(num_domains==1 || mp.only3D_domain_nr)
+      {
+          ret.SetSize(1);
+          // no need to divide mesh, just fill in meshing data
+          ret[0].domain = 1;
+          if(mp.only3D_domain_nr)
+              ret[0].domain = mp.only3D_domain_nr;
+
+          ret[0].mesh.reset(&mesh); // careful, this unique_ptr must not delete &mesh! (it will be released in MergeMeshes after meshing)
+          ret[0].mp = mp;
+          return ret;
+      }
+      ret.SetSize(num_domains);
+
+      Array<Array<PointIndex, PointIndex>> ipmap;
+      ipmap.SetSize(num_domains);
+      auto dim = mesh.GetDimension();
+      auto num_points = mesh.GetNP();
+      auto num_facedescriptors = mesh.GetNFD();
+
+      for(auto i : Range(ret))
+      {
+          auto & md = ret[i];
+          md.domain = i+1;
+
+          md.mp = mp;
+          md.mp.maxh = min2 (mp.maxh, mesh.MaxHDomain(md.domain));
+
+          ret[i].mesh = make_unique<Mesh>();
+          auto & m = *ret[i].mesh;
+
+          m.SetLocalH(mesh.GetLocalH());
+
+          ipmap[i].SetSize(num_points);
+          ipmap[i] = PointIndex::INVALID;
+          m.SetDimension( mesh.GetDimension() );
+          m.SetGeometry( mesh.GetGeometry() );
+
+          for(auto i : Range(1, num_facedescriptors+1))
+              m.AddFaceDescriptor( mesh.GetFaceDescriptor(i) );
+      }
+
+      // mark used points for each domain, add surface elements (with wrong point numbers) to domain mesh
+      for(const auto & sel : mesh.SurfaceElements())
+      {
+        const auto & fd = mesh.GetFaceDescriptor(sel.GetIndex());
+        int dom_in  = fd.DomainIn();
+        int dom_out = fd.DomainOut();
+
+        for( auto dom : {dom_in, dom_out} )
+        {
+            if(dom==0)
+              continue;
+
+            auto & sels = ret[dom-1].mesh->SurfaceElements();
+            for(auto pi : sel.PNums())
+                ipmap[dom-1][pi] = 1;
+            sels.Append(sel);
+        }
+      }
+
+      // mark locked/fixed points for each domain TODO: domain bounding box to add only relevant points?
+      for(auto pi : mesh.LockedPoints())
+        for(auto i : Range(ret))
+          ipmap[i][pi] = 1;
+
+      // add used points to domain mesh, build point mapping
+      for(auto i : Range(ret))
+      {
+          auto & m = *ret[i].mesh;
+          auto & pmap = ret[i].pmap;
+
+          for(auto pi : Range(ipmap[i]))
+            if(ipmap[i][pi])
+            {
+              auto pi_new = m.AddPoint( mesh[pi] );
+              ipmap[i][pi] = pi_new;
+              pmap.Append( pi );
+            }
+      }
+
+      // add segmetns
+      for(auto i : Range(ret))
+      {
+          auto & imap = ipmap[i];
+          auto & m = *ret[i].mesh;
+          for(auto seg : mesh.LineSegments())
+            if(imap[seg[0]].IsValid() && imap[seg[1]].IsValid())
+              {
+                  seg[0] = imap[seg[0]];
+                  seg[1] = imap[seg[1]];
+                  m.AddSegment(seg);
+              }
+      }
+
+      auto & identifications = mesh.GetIdentifications();
+
+      for(auto i : Range(ret))
+      {
+          auto & m = *ret[i].mesh;
+          auto & imap = ipmap[i];
+          auto nmax = identifications.GetMaxNr ();
+          auto & m_ident = m.GetIdentifications();
+
+          for (auto & sel : m.SurfaceElements())
+            for(auto & pi : sel.PNums())
+              pi = imap[pi];
+
+          for(auto n : Range(1,nmax+1))
+          {
+              NgArray<INDEX_2> pairs;
+              identifications.GetPairs(n, pairs);
+
+              for(auto pair : pairs)
+              {
+                  auto pi0 = imap[pair[0]];
+                  auto pi1 = imap[pair[1]];
+                  if(!pi0.IsValid() || !pi1.IsValid())
+                      continue;
+
+                  if(pi1<pi0)
+                      Swap(pi0,pi1);
+                  m_ident.Add(pi0, pi1, n);
+              }
+              m_ident.SetType( n, identifications.GetType(n) );
+          }
+      }
+      return ret;
+  }
+
+  void CloseOpenQuads( MeshingData & md)
+  {
+    auto & mesh = *md.mesh;
+    auto domain = md.domain;
+    MeshingParameters & mp = md.mp;
 
-     MeshingParameters mp = c_mp; // copy mp to change them here
     int oldne;
-     int meshed;
-
-     NgArray<INDEX_2> connectednodes;
-     
-     if (!mesh3d.HasLocalHFunction()) mesh3d.CalcLocalH(mp.grading);
-
-     mesh3d.Compress();
-
-     //  mesh3d.PrintMemInfo (cout);
-
-     if (mp.checkoverlappingboundary)
-        if (mesh3d.CheckOverlappingBoundary())
-           throw NgException ("Stop meshing since boundary mesh is overlapping");
-
-     int nonconsist = 0;
-     for (int k = 1; k <= mesh3d.GetNDomains(); k++)
-     {
-       if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
-	 continue;
-        PrintMessage (3, "Check subdomain ", k, " / ", mesh3d.GetNDomains());
-
-        mesh3d.FindOpenElements(k);
-
-        /*
-        bool res = mesh3d.CheckOverlappingBoundary();
-        if (res)
-        {
-        PrintError ("Surface is overlapping !!");
-        nonconsist = 1;
-        }
-        */
-
-        bool res = (mesh3d.CheckConsistentBoundary() != 0);
-        if (res)
-        {
-           PrintError ("Surface mesh not consistent");
-           nonconsist = 1;
-        }
-     }
-
-     if (nonconsist)
-     {
-        PrintError ("Stop meshing since surface mesh not consistent");
-        throw NgException ("Stop meshing since surface mesh not consistent");
-     }
-
-     double globmaxh = mp.maxh;
-
-     for (int k = 1; k <= mesh3d.GetNDomains(); k++)
-       {
-	 if(mp.only3D_domain_nr && mp.only3D_domain_nr !=k)
-	   continue;
     if (multithread.terminate)
-           break;
+      return;
     
-	 PrintMessage (2, "");
+    mesh.CalcSurfacesOfNode();
-	 PrintMessage (1, "Meshing subdomain ", k, " of ", mesh3d.GetNDomains());
+    mesh.FindOpenElements(domain);
-	 (*testout) << "Meshing subdomain " << k << endl;
-	 
-	 mp.maxh = min2 (globmaxh, mesh3d.MaxHDomain(k));
-
-	 mesh3d.CalcSurfacesOfNode();
-	 mesh3d.FindOpenElements(k);
-	 
-	 if (!mesh3d.GetNOpenElements())
-           continue;
-	 
-	 
-
-	 Box<3> domain_bbox( Box<3>::EMPTY_BOX ); 
-	 
-	 for (SurfaceElementIndex sei = 0; sei < mesh3d.GetNSE(); sei++)
-	   {
-	     const Element2d & el = mesh3d[sei];
-	     if (el.IsDeleted() ) continue;
-
-	     if (mesh3d.GetFaceDescriptor(el.GetIndex()).DomainIn() == k ||
-		 mesh3d.GetFaceDescriptor(el.GetIndex()).DomainOut() == k)
-	       
-	       for (int j = 0; j < el.GetNP(); j++)
-		 domain_bbox.Add (mesh3d[el[j]]);
-	   }
-	 domain_bbox.Increase (0.01 * domain_bbox.Diam());
     
+    if (!mesh.GetNOpenElements())
+      return;
 
     for (int qstep = 0; qstep <= 3; qstep++)
-           // for (int qstep = 0; qstep <= 0; qstep++)  // for hex-filling
      {
        if (qstep == 0 && !mp.try_hexes) continue;
        
-	    // cout << "openquads = " << mesh3d.HasOpenQuads() << endl;
+       if (mesh.HasOpenQuads())
-	    if (mesh3d.HasOpenQuads())
          {
            string rulefile = ngdir;
            
@@ -117,24 +194,19 @@ namespace netgen
            switch (qstep)
              {
              case 0:
-                    rulefile = "/Users/joachim/gitlab/netgen/rules/hexa.rls";
                rulep = hexrules;
                break;
              case 1:
-		    rulefile += "/rules/prisms2.rls";
                rulep = prismrules2;
                break;
              case 2: // connect pyramid to triangle
-		    rulefile += "/rules/pyramids2.rls";
                rulep = pyramidrules2;
                break;
              case 3: // connect to vis-a-vis point
-		    rulefile += "/rules/pyramids.rls";
                rulep = pyramidrules;
                break;
              }
            
-                // Meshing3 meshing(rulefile);
            Meshing3 meshing(rulep); 
            
            MeshingParameters mpquad = mp;
@@ -145,161 +217,172 @@ namespace netgen
            mpquad.check_impossible = qstep == 1;   // for prisms only (air domain in trafo)
            
            
-		// for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
+           for (PointIndex pi : mesh.Points().Range())
-                for (PointIndex pi : mesh3d.Points().Range())
+             meshing.AddPoint (mesh[pi], pi);
-		  meshing.AddPoint (mesh3d[pi], pi);
 
-                /*
+           NgArray<INDEX_2> connectednodes;
-		mesh3d.GetIdentifications().GetPairs (0, connectednodes);
+           for (int nr = 1; nr <= mesh.GetIdentifications().GetMaxNr(); nr++)
-		for (int i = 1; i <= connectednodes.Size(); i++)
+             if (mesh.GetIdentifications().GetType(nr) != Identifications::PERIODIC)
-		  meshing.AddConnectedPair (connectednodes.Get(i));
-                */
-                for (int nr = 1; nr <= mesh3d.GetIdentifications().GetMaxNr(); nr++)
-                  if (mesh3d.GetIdentifications().GetType(nr) != Identifications::PERIODIC)
                {
-                      mesh3d.GetIdentifications().GetPairs (nr, connectednodes);
+                 mesh.GetIdentifications().GetPairs (nr, connectednodes);
                  for (auto pair : connectednodes)
                    meshing.AddConnectedPair (pair);
                }
+           // for (auto pair : md.connected_pairs)
+           //     meshing.AddConnectedPair (pair);
            
-		for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
+           for (int i = 1; i <= mesh.GetNOpenElements(); i++)
              {
-		    Element2d hel = mesh3d.OpenElement(i);
+               Element2d hel = mesh.OpenElement(i);
                meshing.AddBoundaryElement (hel);
              }
            
-		oldne = mesh3d.GetNE();
+           oldne = mesh.GetNE();
            
-		meshing.GenerateMesh (mesh3d, mpquad);
+           meshing.GenerateMesh (mesh, mpquad);
            
-		for (int i = oldne + 1; i <= mesh3d.GetNE(); i++)
+           for (int i = oldne + 1; i <= mesh.GetNE(); i++)
-		  mesh3d.VolumeElement(i).SetIndex (k);
+             mesh.VolumeElement(i).SetIndex (domain);
            
            (*testout) 
-		  << "mesh has " << mesh3d.GetNE() << " prism/pyramid elements" << endl;
+             << "mesh has " << mesh.GetNE() << " prism/pyramid elements" << endl;
            
-		mesh3d.FindOpenElements(k);
+           mesh.FindOpenElements(domain);
          }
      }
    
 
-        if (mesh3d.HasOpenQuads())
+   if (mesh.HasOpenQuads())
    {
       PrintSysError ("mesh has still open quads");
       throw NgException ("Stop meshing since too many attempts");
       // return MESHING3_GIVEUP;
    }
-
-
-        if (mp.delaunay && mesh3d.GetNOpenElements())
-        {
-           Meshing3 meshing((const char**)NULL);
-
-           mesh3d.FindOpenElements(k);
-
-           /*
-           for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
-              meshing.AddPoint (mesh3d[pi], pi);
-           */
-           for (PointIndex pi : mesh3d.Points().Range())
-              meshing.AddPoint (mesh3d[pi], pi);
-
-           for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
-              meshing.AddBoundaryElement (mesh3d.OpenElement(i));
-
-           oldne = mesh3d.GetNE();
-
-           meshing.Delaunay (mesh3d, k, mp);
-
-           for (int i = oldne + 1; i <= mesh3d.GetNE(); i++)
-              mesh3d.VolumeElement(i).SetIndex (k);
-
-           PrintMessage (3, mesh3d.GetNP(), " points, ",
-              mesh3d.GetNE(), " elements");
   }
 
 
+  void MeshDomain( MeshingData & md)
+  {
+    auto & mesh = *md.mesh;
+    auto domain = md.domain;
+    MeshingParameters & mp = md.mp;
+
+    mesh.CalcSurfacesOfNode();
+    mesh.FindOpenElements(md.domain);
+
+    md.meshing = make_unique<Meshing3>(nullptr);
+    for (PointIndex pi : mesh.Points().Range())
+       md.meshing->AddPoint (mesh[pi], pi);
+
+    for (int i = 1; i <= mesh.GetNOpenElements(); i++)
+       md.meshing->AddBoundaryElement (mesh.OpenElement(i));
+
+   if (mp.delaunay && mesh.GetNOpenElements())
+   {
+      int oldne = mesh.GetNE();
+
+      md.meshing->Delaunay (mesh, domain, mp);
+
+      for (int i = oldne + 1; i <= mesh.GetNE(); i++)
+         mesh.VolumeElement(i).SetIndex (domain);
+
+      PrintMessage (3, mesh.GetNP(), " points, ",
+         mesh.GetNE(), " elements");
+   }
+
+   Box<3> domain_bbox( Box<3>::EMPTY_BOX ); 
+   
+   for (auto & sel : mesh.SurfaceElements())
+     {
+       if (sel.IsDeleted() ) continue;
+
+       for (auto pi : sel.PNums())
+         domain_bbox.Add (mesh[pi]);
+     }
+   domain_bbox.Increase (0.01 * domain_bbox.Diam());
+
+   mesh.FindOpenElements(domain);
+
    int cntsteps = 0;
-        if (mesh3d.GetNOpenElements())
+   int meshed;
+   if (mesh.GetNOpenElements())
      do
        {
          if (multithread.terminate)
            break;
          
-              mesh3d.FindOpenElements(k);
+         mesh.FindOpenElements(domain);
-              PrintMessage (5, mesh3d.GetNOpenElements(), " open faces");
+         PrintMessage (5, mesh.GetNOpenElements(), " open faces");
+         GetOpenElements( mesh, domain )->Save("open_"+ToString(cntsteps)+".vol");
          cntsteps++;
 
+
          if (cntsteps > mp.maxoutersteps) 
             throw NgException ("Stop meshing since too many attempts");
 
-              string rulefile = ngdir + "/tetra.rls";
          PrintMessage (1, "start tetmeshing");
 
-              //	  Meshing3 meshing(rulefile);
          Meshing3 meshing(tetrules);
 
-              NgArray<int, PointIndex::BASE> glob2loc(mesh3d.GetNP());
+         Array<PointIndex, PointIndex> glob2loc(mesh.GetNP());
-              glob2loc = -1;
+         glob2loc = PointIndex::INVALID;
 
-              // for (PointIndex pi = mesh3d.Points().Begin(); pi < mesh3d.Points().End(); pi++)
+         for (PointIndex pi : mesh.Points().Range())
-              for (PointIndex pi : mesh3d.Points().Range())
+           if (domain_bbox.IsIn (mesh[pi]))
-                if (domain_bbox.IsIn (mesh3d[pi]))
+               glob2loc[pi] = meshing.AddPoint (mesh[pi], pi);
-                    glob2loc[pi] = 
-                    meshing.AddPoint (mesh3d[pi], pi);
 
-              for (int i = 1; i <= mesh3d.GetNOpenElements(); i++)
+         for (auto sel : mesh.OpenElements() )
          {
-                 Element2d hel = mesh3d.OpenElement(i);
+           for(auto & pi : sel.PNums())
-                 for (int j = 0; j < hel.GetNP(); j++)
+               pi = glob2loc[pi];
-                    hel[j] = glob2loc[hel[j]];
+           meshing.AddBoundaryElement (sel);
-                 meshing.AddBoundaryElement (hel);
-                 // meshing.AddBoundaryElement (mesh3d.OpenElement(i));
          }
 
-              oldne = mesh3d.GetNE();
+         int oldne = mesh.GetNE();
 
          mp.giveuptol = 15 + 10 * cntsteps; 
          mp.sloppy = 5;
-              meshing.GenerateMesh (mesh3d, mp);
+         meshing.GenerateMesh (mesh, mp);
          
-              for (ElementIndex ei = oldne; ei < mesh3d.GetNE(); ei++)
+         for (ElementIndex ei = oldne; ei < mesh.GetNE(); ei++)
-                 mesh3d[ei].SetIndex (k);
+            mesh[ei].SetIndex (domain);
          
 
-              mesh3d.CalcSurfacesOfNode();
+         mesh.CalcSurfacesOfNode();
-              mesh3d.FindOpenElements(k);
+         mesh.FindOpenElements(domain);
 
          // teterrpow = 2;
-              if (mesh3d.GetNOpenElements() != 0)
+         if (mesh.GetNOpenElements() != 0)
          {
             meshed = 0;
-                 PrintMessage (5, mesh3d.GetNOpenElements(), " open faces found");
+            PrintMessage (5, mesh.GetNOpenElements(), " open faces found");
 
             MeshOptimize3d optmesh(mp);
 
             const char * optstr = "mcmstmcmstmcmstmcm";
             for (size_t j = 1; j <= strlen(optstr); j++)
             {
-                    mesh3d.CalcSurfacesOfNode();
+               mesh.CalcSurfacesOfNode();
-                    mesh3d.FreeOpenElementsEnvironment(2);
+               mesh.FreeOpenElementsEnvironment(2);
-                    mesh3d.CalcSurfacesOfNode();
+               mesh.CalcSurfacesOfNode();
 
                switch (optstr[j-1])
                {
-                    case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
+               case 'c': optmesh.CombineImprove(mesh, OPT_REST); break;
-                    case 'd': optmesh.SplitImprove(mesh3d, OPT_REST); break;
+               case 'd': optmesh.SplitImprove(mesh, OPT_REST); break;
-                    case 's': optmesh.SwapImprove(mesh3d, OPT_REST); break;
+               case 's': optmesh.SwapImprove(mesh, OPT_REST); break;
-                    case 't': optmesh.SwapImprove2(mesh3d, OPT_REST); break;
+               case 't': optmesh.SwapImprove2(mesh, OPT_REST); break;
-                    case 'm': mesh3d.ImproveMesh(mp, OPT_REST); break;
+               case 'm': mesh.ImproveMesh(mp, OPT_REST); break;
                }	  
 
             }
 
-                 mesh3d.FindOpenElements(k);	      
+            mesh.FindOpenElements();	      
             PrintMessage (3, "Call remove problem");
-                 RemoveProblem (mesh3d, k);
+            // mesh.Save("before_remove.vol");
-                 mesh3d.FindOpenElements(k);
+            RemoveProblem (mesh, domain);
+            // mesh.Save("after_remove.vol");
+            mesh.FindOpenElements();
          }
          else
            {
@@ -309,11 +392,106 @@ namespace netgen
        }
      while (!meshed);
    
-        PrintMessage (1, mesh3d.GetNP(), " points, ",
+     {
-                      mesh3d.GetNE(), " elements");
+        PrintMessage (3, "Check subdomain ", domain, " / ", mesh.GetNDomains());
+
+        mesh.FindOpenElements(domain);
+
+        bool res = (mesh.CheckConsistentBoundary() != 0);
+        if (res)
+        {
+           PrintError ("Surface mesh not consistent");
+           throw NgException ("Stop meshing since surface mesh not consistent");
+        }
+     }
   }
 
-     mp.maxh = globmaxh;
+  void MergeMeshes( Mesh & mesh, Array<MeshingData> & md )
+  {
+     // todo: optimize: count elements, alloc all memory, copy vol elements in parallel
+     static Timer t("MergeMeshes"); RegionTimer rt(t);
+     if(md.Size()==1)
+     {
+         // assume that mesh was never divided, no need to do anything
+         if(&mesh != md[0].mesh.get())
+             throw Exception("Illegal Mesh pointer in MeshingData");
+
+         md[0].mesh.release();
+         return;
+     }
+
+     for(auto & m_ : md)
+     {
+         auto first_new_pi = m_.pmap.Range().Next();
+         auto & m = *m_.mesh;
+         Array<PointIndex, PointIndex> pmap(m.Points().Size());
+         for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
+             pmap[pi] = m_.pmap[pi];
+
+         for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
+             pmap[pi] = mesh.AddPoint(m[pi]);
+
+
+         for ( auto el : m.VolumeElements() )
+         {
+             for (auto i : Range(el.GetNP()))
+                 el[i] = pmap[el[i]];
+             el.SetIndex(m_.domain);
+             mesh.AddVolumeElement(el);
+         }
+     }
+  }
+
+  void MergeMeshes( Mesh & mesh, FlatArray<Mesh> meshes, PointIndex first_new_pi )
+  {
+     // todo: optimize: count elements, alloc all memory, copy vol elements in parallel
+     static Timer t("MergeMeshes"); RegionTimer rt(t);
+     for(auto & m : meshes)
+     {
+         Array<PointIndex, PointIndex> pmap(m.Points().Size());
+         for(auto pi : Range(PointIndex(PointIndex::BASE), first_new_pi))
+             pmap[pi] = pi;
+
+         for (auto pi : Range(first_new_pi, m.Points().Range().Next()))
+             pmap[pi] = mesh.AddPoint(m[pi]);
+
+
+         for ( auto el : m.VolumeElements() )
+         {
+             for (auto i : Range(el.GetNP()))
+                 el[i] = pmap[el[i]];
+             mesh.AddVolumeElement(el);
+         }
+     }
+  }
+
+  // extern double teterrpow; 
+  MESHING3_RESULT MeshVolume (const MeshingParameters & mp, Mesh& mesh3d)
+  {
+    static Timer t("MeshVolume"); RegionTimer reg(t);
+
+     mesh3d.Compress();
+
+     if (mp.checkoverlappingboundary)
+        if (mesh3d.CheckOverlappingBoundary())
+           throw NgException ("Stop meshing since boundary mesh is overlapping");
+
+
+     if(mesh3d.GetNDomains()==0)
+         return MESHING3_OK;
+
+     if (!mesh3d.HasLocalHFunction())
+         mesh3d.CalcLocalH(mp.grading);
+
+     auto md = DivideMesh(mesh3d, mp);
+
+     ParallelFor( md.Range(), [&](int i)
+       {
+         CloseOpenQuads( md[i] );
+         MeshDomain(md[i]);
+       });
+
+     MergeMeshes(mesh3d, md);
 
      MeshQuality3d (mesh3d);
 
@@ -321,325 +499,6 @@ namespace netgen
   }  
 
 
-
-
-  /*
-
-
-  MESHING3_RESULT MeshVolumeOld (MeshingParameters & mp, Mesh& mesh3d)
-  {
-  int i, k, oldne;
-
-
-  int meshed;
-  int cntsteps; 
-
-
-  PlotStatistics3d * pstat;
-  if (globflags.GetNumFlag("silentflag", 1) <= 2)
-  pstat = new XPlotStatistics3d;
-  else
-  pstat = new TerminalPlotStatistics3d;
-
-  cntsteps = 0;
-  do
-  {
-  cntsteps++;
-  if (cntsteps > mp.maxoutersteps) 
-  {
-  return MESHING3_OUTERSTEPSEXCEEDED;
-  }
-
-
-  int noldp = mesh3d.GetNP();
-      
-      
-  if ( (cntsteps == 1) && globflags.GetDefineFlag ("delaunay"))
-  {
-  cntsteps ++;
-
-  mesh3d.CalcSurfacesOfNode();
-
-
-  for (k = 1; k <= mesh3d.GetNDomains(); k++)
-  {
-  Meshing3 meshing(NULL, pstat);
-
-  mesh3d.FindOpenElements(k);
-	      
-  for (i = 1; i <= noldp; i++)
-  meshing.AddPoint (mesh3d.Point(i), i);
-	      
-  for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
-  {
-  if (mesh3d.OpenElement(i).GetIndex() == k)
-  meshing.AddBoundaryElement (mesh3d.OpenElement(i));
-  }
-	      
-  oldne = mesh3d.GetNE();
-  if (globflags.GetDefineFlag ("blockfill"))
-  {
-  if (!globflags.GetDefineFlag ("localh"))
-  meshing.BlockFill 
-  (mesh3d, mp.h * globflags.GetNumFlag ("relblockfillh", 1));
-  else
-  meshing.BlockFillLocalH (mesh3d);
-  }
-	      
-  MeshingParameters mpd;
-  meshing.Delaunay (mesh3d, mpd);
-
-  for (i = oldne + 1; i <= mesh3d.GetNE(); i++)
-  mesh3d.VolumeElement(i).SetIndex (k);
-  }
-  }
-
-  noldp = mesh3d.GetNP();
-
-  mesh3d.CalcSurfacesOfNode();
-  mesh3d.FindOpenElements();
-  for (k = 1; k <= mesh3d.GetNDomains(); k++)
-  {
-  Meshing3 meshing(globflags.GetStringFlag ("rules3d", NULL), pstat);
-      
-  Point3d pmin, pmax;
-  mesh3d.GetBox (pmin, pmax, k);
-	  
-  rot.SetCenter (Center (pmin, pmax));
-
-  for (i = 1; i <= noldp; i++)
-  meshing.AddPoint (mesh3d.Point(i), i);
-
-  for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
-  {
-  if (mesh3d.OpenElement(i).GetIndex() == k)
-  meshing.AddBoundaryElement (mesh3d.OpenElement(i));
-  }
-
-  oldne = mesh3d.GetNE();
-
-
-  if ( (cntsteps == 1) && globflags.GetDefineFlag ("blockfill"))
-  {
-  if (!globflags.GetDefineFlag ("localh"))
-  {
-  meshing.BlockFill 
-  (mesh3d, 
-  mp.h * globflags.GetNumFlag ("relblockfillh", 1));
-  }
-  else
-  {
-  meshing.BlockFillLocalH (mesh3d);
-  }
-  }
-
-
-  mp.giveuptol = int(globflags.GetNumFlag ("giveuptol", 15));
-
-  meshing.GenerateMesh (mesh3d, mp);
-
-  for (i = oldne + 1; i <= mesh3d.GetNE(); i++)
-  mesh3d.VolumeElement(i).SetIndex (k);
-  }
-
-
-
-  mesh3d.CalcSurfacesOfNode();
-  mesh3d.FindOpenElements();
-      
-  teterrpow = 2;
-  if (mesh3d.GetNOpenElements() != 0)
-  {
-  meshed = 0;
-  (*mycout) << "Open elements found, old" << endl;
-  const char * optstr = "mcmcmcmcm";
-  int j;
-  for (j = 1; j <= strlen(optstr); j++)
-  switch (optstr[j-1])
-  {
-  case 'c': mesh3d.CombineImprove(); break;
-  case 'd': mesh3d.SplitImprove(); break;
-  case 's': mesh3d.SwapImprove(); break;
-  case 'm': mesh3d.ImproveMesh(2); break;
-  }	  
-	  
-  (*mycout) << "Call remove" << endl;
-  RemoveProblem (mesh3d);
-  (*mycout) << "Problem removed" << endl;
-  }
-  else
-  meshed = 1;
-  }
-  while (!meshed);
-
-  MeshQuality3d (mesh3d);
-
-  return MESHING3_OK;
-  }  
-
-  */
-
-
-
-
-  /*
-  MESHING3_RESULT MeshMixedVolume(MeshingParameters & mp, Mesh& mesh3d)
-  {
-    int i, j;
-    MESHING3_RESULT res;
-    Point3d pmin, pmax;
-
-    mp.giveuptol = 10;
-    mp.baseelnp = 4;
-    mp.starshapeclass = 100;
-
-    //  TerminalPlotStatistics3d pstat;
-  
-    Meshing3 meshing1("pyramids.rls");
-    for (i = 1; i <= mesh3d.GetNP(); i++)
-      meshing1.AddPoint (mesh3d.Point(i), i);
-
-    mesh3d.FindOpenElements();
-    for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
-      if (mesh3d.OpenElement(i).GetIndex() == 1)
-	meshing1.AddBoundaryElement (mesh3d.OpenElement(i));
-
-    res = meshing1.GenerateMesh (mesh3d, mp);
-
-    mesh3d.GetBox (pmin, pmax);
-    PrintMessage (1, "Mesh pyramids, res = ", res);
-    if (res)
-      exit (1);
-
-
-    for (i = 1; i <= mesh3d.GetNE(); i++)
-      mesh3d.VolumeElement(i).SetIndex (1);
-
-    // do delaunay
-  
-    mp.baseelnp = 0;
-    mp.starshapeclass = 5;
-
-    Meshing3 meshing2(NULL);
-    for (i = 1; i <= mesh3d.GetNP(); i++)
-      meshing2.AddPoint (mesh3d.Point(i), i);
-    
-    mesh3d.FindOpenElements();
-    for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
-      if (mesh3d.OpenElement(i).GetIndex() == 1)
-	meshing2.AddBoundaryElement (mesh3d.OpenElement(i));
-
-    MeshingParameters mpd;
-    meshing2.Delaunay (mesh3d, mpd);
-
-    for (i = 1; i <= mesh3d.GetNE(); i++)
-      mesh3d.VolumeElement(i).SetIndex (1);
-
-
-    mp.baseelnp = 0;
-    mp.giveuptol = 10;
-
-    for (int trials = 1; trials <= 50; trials++)
-      {
-	if (multithread.terminate)
-	  return MESHING3_TERMINATE;
-
-	Meshing3 meshing3("tetra.rls");
-	for (i = 1; i <= mesh3d.GetNP(); i++)
-	  meshing3.AddPoint (mesh3d.Point(i), i);
-      
-	mesh3d.FindOpenElements();
-	for (i = 1; i <= mesh3d.GetNOpenElements(); i++)
-	  if (mesh3d.OpenElement(i).GetIndex() == 1)
-	    meshing3.AddBoundaryElement (mesh3d.OpenElement(i));
-      
-	if (trials > 1)
-	  CheckSurfaceMesh2 (mesh3d);
-	res = meshing3.GenerateMesh (mesh3d, mp);
-      
-	for (i = 1; i <= mesh3d.GetNE(); i++)
-	  mesh3d.VolumeElement(i).SetIndex (1);
-
-	if (res == 0) break;
-
-
-
-	for (i = 1; i <= mesh3d.GetNE(); i++)
-	  {
-	    const Element & el = mesh3d.VolumeElement(i);
-	    if (el.GetNP() != 4)
-	      {
-		for (j = 1; j <= el.GetNP(); j++)
-		  mesh3d.AddLockedPoint (el.PNum(j));
-	      }
-	  }
-
-	mesh3d.CalcSurfacesOfNode();
-	mesh3d.FindOpenElements();
-
-	MeshOptimize3d optmesh;
-
-	teterrpow = 2;
-	const char * optstr = "mcmcmcmcm";
-	for (j = 1; j <= strlen(optstr); j++)
-	  switch (optstr[j-1])
-	    {
-	    case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
-	    case 'd': optmesh.SplitImprove(mesh3d); break;
-	    case 's': optmesh.SwapImprove(mesh3d); break;
-	    case 'm': mesh3d.ImproveMesh(); break;
-	    }	  
-	        
-	RemoveProblem (mesh3d);
-      }
-
-
-    PrintMessage (1, "Meshing tets, res = ", res);
-    if (res)
-      {
-	mesh3d.FindOpenElements();
-	PrintSysError (1, "Open elements: ", mesh3d.GetNOpenElements());
-	exit (1);
-      }
-
-
-  
-    for (i = 1; i <= mesh3d.GetNE(); i++)
-      {
-	const Element & el = mesh3d.VolumeElement(i);
-	if (el.GetNP() != 4)
-	  {
-	    for (j = 1; j <= el.GetNP(); j++)
-	      mesh3d.AddLockedPoint (el.PNum(j));
-	  }
-      }
-  
-    mesh3d.CalcSurfacesOfNode();
-    mesh3d.FindOpenElements();
-  
-    MeshOptimize3d optmesh;
-
-    teterrpow = 2;
-    const char * optstr = "mcmcmcmcm";
-    for (j = 1; j <= strlen(optstr); j++)
-      switch (optstr[j-1])
-	{
-	case 'c': optmesh.CombineImprove(mesh3d, OPT_REST); break;
-	case 'd': optmesh.SplitImprove(mesh3d); break;
-	case 's': optmesh.SwapImprove(mesh3d); break;
-	case 'm': mesh3d.ImproveMesh(); break;
-	}	  
-
-
-    return MESHING3_OK;
-  }
-*/
-
-
-
-
-
-
   MESHING3_RESULT OptimizeVolume (const MeshingParameters & mp, 
 				  Mesh & mesh3d)
     //				  const CSGeometry * geometry)

libsrc/meshing/meshing3.cpp

@@ -1144,11 +1144,18 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
 
   if (mp.maxh < maxh) maxh = mp.maxh;
 
+  auto loch_ptr = mesh.LocalHFunction().Copy(bbox);
+  auto & loch = *loch_ptr;
+
   bool changed;
+  static Timer t1("loop1");
+  t1.Start();
   do 
     {
-      mesh.LocalHFunction().ClearFlags();
+      loch.ClearFlags();
 
+      static Timer tbox("adfront-bbox");
+      tbox.Start();
       for (int i = 1; i <= adfront->GetNF(); i++)
 	{
 	  const MiniElement2d & el = adfront->GetFace(i);
@@ -1161,26 +1168,28 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
 	  double filld = filldist * bbox.Diam();
 	  bbox.Increase (filld);
       
-      	  mesh.LocalHFunction().CutBoundary (bbox); // .PMin(), bbox.PMax());
+      	  loch.CutBoundary (bbox); // .PMin(), bbox.PMax());
 	}
+      tbox.Stop();
 
       //      locadfront = adfront;
-      mesh.LocalHFunction().FindInnerBoxes (adfront, NULL);
+      loch.FindInnerBoxes (adfront, NULL);
 
       npoints.SetSize(0);
-      mesh.LocalHFunction().GetInnerPoints (npoints);
+      loch.GetInnerPoints (npoints);
 
       changed = false;
       for (int i = 1; i <= npoints.Size(); i++)
 	{
-	  if (mesh.LocalHFunction().GetH(npoints.Get(i)) > 1.5 * maxh)
+	  if (loch.GetH(npoints.Get(i)) > 1.5 * maxh)
 	    {
-	      mesh.LocalHFunction().SetH (npoints.Get(i), maxh);
+	      loch.SetH (npoints.Get(i), maxh);
 	      changed = true;
 	    }
 	}
     }
   while (changed);
+  t1.Stop();
 
   if (debugparam.slowchecks)
     (*testout) << "Blockfill with points: " << endl;
@@ -1208,6 +1217,8 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
 
   // find outer points
   
+  static Timer tloch2("build loch2");
+  tloch2.Start();
   loch2.ClearFlags();
 
   for (int i = 1; i <= adfront->GetNF(); i++)
@@ -1245,6 +1256,7 @@ void Meshing3 :: BlockFillLocalH (Mesh & mesh,
       // loch2.CutBoundary (pmin, pmax);
       loch2.CutBoundary (Box<3> (pmin, pmax)); // pmin, pmax);
     }
+  tloch2.Stop();
 
   // locadfront = adfront;
   loch2.FindInnerBoxes (adfront, NULL);

libsrc/meshing/meshtype.hpp

@@ -852,13 +852,24 @@ namespace netgen
         { _np = np; _typ = typ; _curved = is_curved; }
       // ar & _np & _typ & index & _curved;
       ar.DoPacked (_np, _typ, index, _curved);                
-      if (ar.Input())
-        { np = _np; typ = ELEMENT_TYPE(_typ); is_curved = _curved; }
 
-      /*
+      if (ar.Input())
-      for (size_t i = 0; i < np; i++)
+        {
-        ar & pnum[i];
+          np = _np;
-      */
+          typ = ELEMENT_TYPE(_typ);
+          is_curved = _curved;
+          flags.marked = 1;
+          flags.badel = 0;
+          flags.reverse = 0;
+          flags.illegal = 0;
+          flags.illegal_valid = 0;
+          flags.badness_valid = 0;
+          flags.refflag = 1;
+          flags.strongrefflag = false;
+          flags.deleted = 0;
+          flags.fixed = 0;
+        }
+
       static_assert(sizeof(int) == sizeof (PointIndex));
       ar.Do( (int*)&pnum[0], np);
     }

tests/pytest/compare_results.py

@@ -15,7 +15,6 @@ def readData(a, files):
     file=[]
     for f in files:
         for t in a[f]:
-            file.append(f)
             if t['ne1d']>0:
                 ne1d.append(t['ne1d'])
             if t['ne2d']>0:
@@ -24,6 +23,7 @@ def readData(a, files):
                 ne3d.append(t['ne3d'])
             if t['total_badness']>0.0:
                 bad.append(t['total_badness'])
+                file.append(f)
                 if 'angles_tet' in t:
                     amin.append(t['angles_tet'][0])
                     amax.append(t['angles_tet'][1])

tests/pytest/test_pickling.py

@@ -87,9 +87,13 @@ def test_pickle_geom2d():
 
 def test_pickle_mesh():
     import netgen.csg as csg
-    geo = csg.CSGeometry()
+    geo1 = csg.CSGeometry()
+    geo2 = csg.CSGeometry()
     brick = csg.OrthoBrick(csg.Pnt(-3,-3,-3), csg.Pnt(3,3,3))
-    mesh = geo.GenerateMesh(maxh=0.2)
+    geo2.Add(brick)
+
+    for geo in [geo1, geo2]:
+        mesh = geo.GenerateMesh(maxh=2)
         assert geo == mesh.GetGeometry()
         dump = pickle.dumps([geo,mesh])
         geo2, mesh2 = pickle.loads(dump)