Merge branch 'occ_zrefinement' into 'master'

ZRefinement for arbitrary geometry class

See merge request jschoeberl/netgen!479

libsrc/general/hashtabl.hpp

@@ -150,6 +150,13 @@ public:
     return cont.Get (bnr, pos);
   }
 
+  T & Get (const INDEX_2 & ahash)
+  {
+    int bnr = HashValue (ahash);
+    int pos = Position (bnr, ahash);
+    return cont.Get (bnr, pos);
+  }
+
   ///
   bool Used (const INDEX_2 & ahash) const
   {
@@ -214,9 +221,14 @@ public:
     int BagNr() const { return bagnr; }
     int Pos() const { return pos; }
 
-    void operator++ (int)
+    Iterator operator++ (int)
+    {
+      Iterator it(ht, bagnr, pos);
+      ++(*this);
+      return it;
+    }
+    Iterator& operator++()
     {
-      // cout << "begin Operator ++: bagnr = " << bagnr << " -  pos = " << pos << endl;
       pos++;
       while (bagnr < ht.GetNBags() && 
 	     pos == ht.GetBagSize(bagnr+1))
@@ -224,7 +236,12 @@ public:
 	  pos = 0;
 	  bagnr++;
 	}
-      // cout << "end Operator ++: bagnr = " << bagnr << " - pos = " << pos << endl;
+      return *this;
+    }
+
+    std::pair<INDEX_2, T> operator*()
+    {
+      return std::make_pair(ht.hash[bagnr][pos], ht.cont[bagnr][pos]);
     }
 
     bool operator != (int i) const
@@ -246,6 +263,18 @@ public:
     return GetNBags();
   }
 
+  Iterator begin () const
+  {
+    Iterator it(*this, 0, -1);
+    it++;
+    return it;
+  }
+
+  int end() const
+  {
+    return GetNBags();
+  }
+
   void GetData (const Iterator & it,
 		INDEX_2 & ahash, T & acont) const
   {

libsrc/general/table.hpp

@@ -201,6 +201,8 @@ public:
   inline const T & Get (int i, int nr) const
     { return ((T*)data.Get(i).col)[nr-1]; }
 
+  inline T & Get (int i, int nr)
+    { return ((T*)data.Get(i).col)[nr-1]; }
 
   /** Returns pointer to the first element in row i. */
   inline const T * GetLine (int i) const

libsrc/meshing/meshclass.cpp

@@ -6213,7 +6213,169 @@ namespace netgen
     */
   }
 
+  void Mesh :: ZRefine(const string& name, const Array<double>& slices)
+  {
+    auto nr = GetIdentifications().GetNr(name);
+    auto& identpts = GetIdentifications().GetIdentifiedPoints();
 
+    UpdateTopology();
+
+    std::map<std::pair<PointIndex, PointIndex>,
+             Array<PointIndex>> inserted_points;
+    BitArray mapped_points(GetNV()+1);
+    mapped_points = false;
+
+    // Add new points
+    for(auto [p1p2, idnr] : identpts)
+      {
+        if(idnr != nr)
+          continue;
+        auto& ipts = inserted_points[{p1p2.I1(), p1p2.I2()}];
+        auto p1 = Point(p1p2.I1());
+        auto p2 = Point(p1p2.I2());
+        ipts.Append(p1p2.I1());
+        mapped_points.SetBit(p1p2.I1());
+        for(auto slice : slices)
+          {
+            auto np = p1 + slice * (p2-p1);
+            auto npi = AddPoint(np);
+            ipts.Append(npi);
+          }
+        ipts.Append(p1p2.I2());
+      }
+
+    // Split segments
+    for(auto si : Range(segments))
+      {
+        auto& seg = segments[si];
+        auto p1 = seg[0];
+        auto p2 = seg[1];
+
+        auto c1 = inserted_points.count({p1, p2});
+        auto c2 = inserted_points.count({p2, p1});
+
+        if(c1 == 0 && c2 == 0)
+          continue;
+
+        if(c2)
+          Swap(p1,p2);
+
+        const auto& ipts = inserted_points[{p1,p2}];
+        if(c2)
+          seg[1] = ipts[ipts.Size()-2];
+        else
+          seg[1] = ipts[1];
+        for(auto i : Range(size_t(1), ipts.Size()-1))
+          {
+            Segment snew = seg;
+            if(c2)
+              {
+                seg[0] = ipts[ipts.Size()-1-i];
+                seg[1] = ipts[ipts.Size()-2-i];
+              }
+            else
+              {
+                snew[0] = ipts[i];
+                snew[1] = ipts[i+1];
+              }
+            AddSegment(snew);
+          }
+      }
+
+    BitArray sel_done(surfelements.Size());
+    sel_done = false;
+
+    // Split surface elements
+    auto p2sel = CreatePoint2SurfaceElementTable();
+    for(const auto& [pair, inserted] : inserted_points)
+      {
+        for(auto si : p2sel[pair.first])
+          {
+            if(sel_done[si])
+              continue;
+            sel_done.SetBit(si);
+            auto sel = surfelements[si];
+            map<PointIndex, Array<PointIndex>> mapped_points;
+            int nmapped = 0;
+            for(auto i : Range(sel.GetNP()))
+              {
+                auto p1 = sel[i];
+                auto p2 = sel[(i+1)%sel.GetNP()];
+                auto c1 = inserted_points.count({p1, p2});
+                auto c2 = inserted_points.count({p2, p1});
+                if(c1 == 0 && c2 == 0)
+                  continue;
+                if(c2)
+                  Swap(p1, p2);
+                auto& ipts = inserted_points[{p1, p2}];
+                auto& a1 = mapped_points[p1];
+                auto& a2 = mapped_points[p2];
+                a1 = ipts.Range(0, ipts.Size()-1);
+                a2 = ipts.Range(1, ipts.Size());
+                nmapped = ipts.Size()-1;
+              }
+            for(auto i : Range(nmapped))
+              {
+                Element2d nsel = sel;
+                for(auto& pi : nsel.PNums())
+                  if(mapped_points.count(pi))
+                    pi = mapped_points[pi][i];
+                AddSurfaceElement(nsel);
+              }
+            if(nmapped)
+              surfelements[si].Delete();
+          }
+      }
+
+    // Split volume elements
+    BitArray vol_done(volelements.Size());
+    vol_done = false;
+    auto p2el = CreatePoint2ElementTable(); // mapped_points);
+    for(const auto& [pair, inserted] : inserted_points)
+      {
+        for(auto ei : p2el[pair.first])
+          {
+            if(vol_done[ei])
+              continue;
+            vol_done.SetBit(ei);
+            auto el = volelements[ei];
+            map<PointIndex, Array<PointIndex>> mapped_points;
+            int nmapped = 0;
+            NgArray<int> eledges;
+            topology.GetElementEdges(ei+1, eledges);
+            for(auto edgei : eledges)
+              {
+                int p1, p2;
+                topology.GetEdgeVertices(edgei, p1, p2);
+                auto c1 = inserted_points.count({p1, p2});
+                auto c2 = inserted_points.count({p2, p1});
+                if(c1 == 0 && c2 == 0)
+                  continue;
+                if(c2)
+                  Swap(p1, p2);
+                auto& ipts = inserted_points[{p1, p2}];
+                auto& a1 = mapped_points[p1];
+                auto& a2 = mapped_points[p2];
+                a1 = ipts.Range(0, ipts.Size()-1);
+                a2 = ipts.Range(1, ipts.Size());
+                nmapped = ipts.Size()-1;
+              }
+
+            for(auto i : Range(nmapped))
+              {
+                Element nel = el;
+                for(auto& pi : nel.PNums())
+                  if(mapped_points.count(pi))
+                    pi = mapped_points[pi][i];
+                AddVolumeElement(nel);
+              }
+            if(nmapped)
+              volelements[ei].Delete();
+          }
+      }
+
+    Compress();
+  }
 
   void Mesh :: RebuildSurfaceElementLists ()
   {

libsrc/meshing/meshclass.hpp

@@ -499,6 +499,7 @@ namespace netgen
     /// Refines mesh and projects points to true surface
     // void Refine (int levels, const CSGeometry * geom);
 
+    void ZRefine(const string& name, const Array<double>& slices);
     
     bool BoundaryEdge (PointIndex pi1, PointIndex pi2) const
     {

libsrc/meshing/python_mesh.cpp

@@ -1127,6 +1127,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
              self.UpdateTopology();
            }),py::call_guard<py::gil_scoped_release>())
 
+    .def("ZRefine", &Mesh::ZRefine)
+
     .def ("SecondOrder", FunctionPointer
           ([](Mesh & self)
            {