fix boundarylayer 2d code (now single line segments, not per face)

libsrc/meshing/boundarylayer.cpp

@@ -436,7 +436,6 @@ namespace netgen
                   }
               }
 
-              // seg.edgenr = topo.GetEdge(segi)+1;
               segments.Append(seg);
           }
       }
@@ -1394,20 +1393,24 @@ namespace netgen
      compress = PointIndex::INVALID;
 
      PointIndex cnt = PointIndex::BASE;
-     for(const auto & seg : mesh.LineSegments())
-        if (seg.si == domain)
-           for (auto pi : {seg[0], seg[1]})
-              if (compress[pi]==PointIndex{PointIndex::INVALID})
-              {
-                 meshing.AddPoint(mesh[pi], pi);
-                 compress[pi] = cnt++;
-              }
 
+     auto p2sel = mesh.CreatePoint2SurfaceElementTable();
+     Array<Segment> domain_segments;
      PointGeomInfo gi;
      gi.trignum = domain;
-     for(const auto & seg : mesh.LineSegments())
+     for(auto seg : mesh.LineSegments())
-        if (seg.si == domain)
+     {
-           meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi);
+         for (auto pi : {seg[0], seg[1]})
+            if (compress[pi]==PointIndex{PointIndex::INVALID})
+            {
+               meshing.AddPoint(mesh[pi], pi);
+               compress[pi] = cnt++;
+            }
+         if(seg.domin == domain)
+             meshing.AddBoundaryElement (compress[seg[0]], compress[seg[1]], gi, gi);
+         if(seg.domout == domain)
+             meshing.AddBoundaryElement (compress[seg[1]], compress[seg[0]], gi, gi);
+     }
 
      auto oldnf = mesh.GetNSE();
      auto res = meshing.GenerateMesh (mesh, mp, mp.maxh, domain);
@@ -1450,17 +1453,20 @@ namespace netgen
      }
 
      mesh.Compress();
+     mesh.CalcSurfacesOfNode();
      mesh.OrderElements();
      mesh.SetNextMajorTimeStamp();
-
   }
 
   int GenerateBoundaryLayer2 (Mesh & mesh, int domain, const Array<double> & thicknesses, bool should_make_new_domain, const Array<int> & boundaries)
   {
+     mesh.GetTopology().SetBuildVertex2Element(true);
+     mesh.UpdateTopology();
+     const auto & line_segments = mesh.LineSegments();
      SegmentIndex first_new_seg = mesh.LineSegments().Range().Next();
 
      int np = mesh.GetNP();
-     int nseg = mesh.GetNSeg();
+     int nseg = line_segments.Size();
      int ne = mesh.GetNSE();
      mesh.UpdateTopology();
 
@@ -1482,23 +1488,10 @@ namespace netgen
 
      auto & meshtopo = mesh.GetTopology();
 
-     Array<SurfaceElementIndex, SegmentIndex> seg2surfel(mesh.GetNSeg());
-     seg2surfel = -1;
-     NgArray<SurfaceElementIndex> temp_els;
-     for(auto si : Range(mesh.LineSegments()))
-     {
-        meshtopo.GetSegmentSurfaceElements ( si+1, temp_els );
-        // NgArray<int> surfeledges;
-        // meshtopo.GetSurfaceElementEdges(si+1, surfeledges);
-        for(auto seli : temp_els)
-           if(mesh[seli].GetIndex() == mesh[si].si)
-              seg2surfel[si] = seli;
-     }
-
      Array<SegmentIndex> segments;
 
     // surface index map
-    Array<int> si_map(mesh.GetNFD()+1);
+    Array<int> si_map(mesh.GetNFD()+2);
     si_map = -1;
 
     int fd_old = mesh.GetNFD();
@@ -1506,12 +1499,14 @@ namespace netgen
     int max_edge_nr = -1;
     int max_domain = -1;
 
-    for(const auto& seg : mesh.LineSegments())
+    for(const auto& seg : line_segments)
     {
       if(seg.epgeominfo[0].edgenr > max_edge_nr)
         max_edge_nr = seg.epgeominfo[0].edgenr;
-      if(seg.si > max_domain)
+      if(seg.domin > max_domain)
-         max_domain = seg.si;
+         max_domain = seg.domin;
+      if(seg.domout > max_domain)
+         max_domain = seg.domout;
     }
 
     int new_domain = max_domain+1;
@@ -1527,95 +1522,43 @@ namespace netgen
        for(auto edgenr : boundaries)
           active_boundaries.SetBit(edgenr);
 
-    for(auto segi : Range(mesh.LineSegments()))
+    for(auto segi : Range(line_segments))
     {
-       const auto seg = mesh[segi];
+       const auto seg = line_segments[segi];
-       if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && seg.si==domain)
+       if(active_boundaries.Test(seg.epgeominfo[0].edgenr) && (seg.domin==domain || seg.domout==domain))
           active_segments.SetBit(segi);
     }
 
-    for(auto segi : Range(mesh.LineSegments()))
     {
-        const auto& seg = mesh[segi];
-        auto si = seg.si;
-
-        if(si_map[si]!=-1)
-           continue;
-
-        if(!active_segments.Test(segi))
-           continue;
-
         FaceDescriptor new_fd(0, 0, 0, -1);
         new_fd.SetBCProperty(new_domain);
         int new_fd_index = mesh.AddFaceDescriptor(new_fd);
-        si_map[si] = new_domain;
         if(should_make_new_domain)
-           mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(si-1));
+           mesh.SetBCName(new_domain-1, "mapped_" + mesh.GetBCName(domain-1));
     }
 
-    for(auto si : Range(mesh.LineSegments()))
+    for(auto segi : Range(line_segments))
       {
-        if(segs_done[si]) continue;
+        if(segs_done[segi]) continue;
-        segs_done.SetBit(si);
+        segs_done.SetBit(segi);
-        const auto& segi = mesh[si];
+        const auto& seg = line_segments[segi];
-        if(si_map[segi.si] == -1) continue;
+        if(seg.domin != domain && seg.domout != domain) continue;
-        if(!active_boundaries.Test(segi.epgeominfo[0].edgenr))
+        if(!active_boundaries.Test(seg.epgeominfo[0].edgenr))
            continue;
-        moved_segs.Append(si);
+        moved_segs.Append(segi);
       }
 
      // calculate growth vectors (average normal vectors of adjacent segments at each point)
      for (auto si : moved_segs)
      {
-       auto & seg = mesh[si];
+       auto & seg = line_segments[si];
-
-       meshtopo.GetSegmentSurfaceElements ( si+1, temp_els );
-       ArrayMem<int, 10> seg_domains;
-
-       temp_els.SetSize(0);
-       if(seg2surfel[si]!=-1)
-          temp_els.Append(seg2surfel[si]);
-
-       int n_temp_els = temp_els.Size();
-       if(n_temp_els==0)
-          continue;
-
-       int dom0 = mesh[temp_els[0]].GetIndex();
-       int dom1 = n_temp_els==2 ? mesh[temp_els[1]].GetIndex() : 0;
-
-       bool in_dom0 = dom0 == domain;
-       bool in_dom1 = dom1 == domain;
-
-       if(!in_dom0 && !in_dom1)
-          continue;
-
-       int side = in_dom0 ? 0 : 1;
-
-       auto & sel = mesh[ temp_els[side] ];
-
-       int domain = sel.GetIndex();
-       Vec<3> pcenter = 0.0;
-       for(auto i : IntRange(sel.GetNP()))
-       {
-          for(auto d : IntRange(3))
-             pcenter[d] += mesh[sel[i]][d];
-       }
-       pcenter = 1.0/sel.GetNP() * pcenter;
 
        auto n = mesh[seg[1]] - mesh[seg[0]];
        n = {-n[1], n[0], 0};
        n.Normalize();
 
-       Vec<3> p0{mesh[seg[0]]};
+       if(seg.domout == domain)
-       Vec<3> p1{mesh[seg[0]]};
+           n = -n;
-
-
-       auto v = pcenter -0.5*(p0+p1);
-
-       const auto Dot = [](Vec<3> a, Vec<3> b)
-       { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; };
-       if(Dot(n, v)<0)
-          n = -1*n;
 
        AddDirection(growthvectors[seg[0]], n);
        AddDirection(growthvectors[seg[1]], n);
@@ -1628,7 +1571,7 @@ namespace netgen
 
      for(auto si : moved_segs)
      {
-        auto current_seg = mesh[si];
+        auto current_seg = line_segments[si];
         auto current_si = si;
 
         auto first = current_seg[0];
@@ -1766,8 +1709,9 @@ namespace netgen
         const auto & seg0 = mesh[segi0];
         const auto & seg1 = mesh[segi1];
 
-        if(seg0.si != seg1.si)
+        if( (seg0.domin != domain && seg0.domout != domain) ||
-           return;
+            (seg1.domin != domain && seg1.domout != domain) )
+            return;
 
         if(segi0 == segi1)
            return;
@@ -1880,7 +1824,7 @@ namespace netgen
      {
         auto p2 = [](Point<3> p) { return Point<2>{p[0], p[1]}; };
 
-        auto seg = mesh[segi];
+        auto seg = line_segments[segi];
         double alpha,beta;
         intersect( p2(mesh[seg[0]]), p2(mesh[seg[0]]+total_thickness*growthvectors[seg[0]]), p2(mesh[seg[1]]), p2(mesh[seg[1]]+total_thickness*growthvectors[seg[1]]), alpha, beta );
 
@@ -1918,13 +1862,13 @@ namespace netgen
      // insert new elements ( and move old ones )
      for(auto si : moved_segs)
      {
-        auto seg = mesh[si];
+        auto seg = line_segments[si];
 
         bool swap = false;
         auto & pm0 = mapto[seg[0]];
         auto & pm1 = mapto[seg[1]];
 
-        auto newindex = si_map[seg.si];
+        auto newindex = si_map[domain];
 
         Segment s = seg;
         s.geominfo[0] = {};
@@ -1939,10 +1883,6 @@ namespace netgen
         s.si = seg.si;
         mesh.AddSegment(s);
 
-        Swap(s[0], s[1]);
-        s.si =  newindex;
-        mesh.AddSegment(s);
-
         for ( auto i : Range(thicknesses))
         {
            PointIndex pi0, pi1, pi2, pi3;
@@ -1980,14 +1920,14 @@ namespace netgen
            newel[1] = pi1;
            newel[2] = pi2;
            newel[3] = pi3;
-           newel.SetIndex(si_map[seg.si]);
+           newel.SetIndex(new_domain);
            newel.GeomInfo() = PointGeomInfo{};
 
-//            if(swap)
+            if(swap)
-//            {
+            {
-//               Swap(newel[0], newel[1]);
+               Swap(newel[0], newel[1]);
-//               Swap(newel[2], newel[3]);
+               Swap(newel[2], newel[3]);
-//            }
+            }
 
            for(auto i : Range(4))
            {
@@ -1998,7 +1938,10 @@ namespace netgen
 
         }
         // segment now adjacent to new 2d-domain!
-        mesh[si].si = si_map[seg.si];
+        if(line_segments[si].domin == domain)
+            line_segments[si].domin = new_domain;
+        if(line_segments[si].domout == domain)
+            line_segments[si].domout = new_domain;
      }
 
      for(auto pi : Range(mapto))
@@ -2043,7 +1986,12 @@ namespace netgen
         }
 
         for(auto segi : moved_segs)
-           mesh[segi].si = domain;
+        {
+            if(mesh[segi].domin == new_domain)
+                mesh[segi].domin = domain;
+            if(mesh[segi].domout == new_domain)
+                mesh[segi].domout = domain;
+        }
 
         mesh.Compress();
         mesh.CalcSurfacesOfNode();

libsrc/meshing/fieldlines.cpp

@@ -77,7 +77,7 @@ namespace netgen
   {
     bool finished = false;
     
-    if(stepcount <= steps)
+    if(stepcount <= steps && stepcount>0)
       {
 	t = startt + c[stepcount-1]*h;
 	val = startval;

libsrc/meshing/python_mesh.cpp

@@ -1191,6 +1191,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m)
     
     .def ("BuildSearchTree", &Mesh::BuildElementSearchTree,py::call_guard<py::gil_scoped_release>())
 
+    .def ("BoundaryLayer2", GenerateBoundaryLayer2, py::arg("domain"), py::arg("thicknesses"), py::arg("make_new_domain")=true, py::arg("boundaries")=Array<int>{})
     .def ("BoundaryLayer", [](Mesh & self, variant<string, int> boundary,
                               variant<double, py::list> thickness,
                               string material,

libsrc/occ/occgeom.hpp

@@ -202,6 +202,12 @@ namespace netgen
     Mesh::GEOM_TYPE GetGeomType() const override
     { return Mesh::GEOM_OCC; }
 
+    void SetDimension(int dim)
+    {
+        dimension = dim;
+        BuildFMap();
+    }
+
     void SetOCCParameters(const OCCParameters& par)
     { occparam = par; }
 

libsrc/occ/python_occ.cpp

@@ -134,7 +134,7 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                   }), py::arg("shape"),
          "Create Netgen OCCGeometry from existing TopoDS_Shape")
     
-    .def(py::init([] (const string& filename)
+    .def(py::init([] (const string& filename, int dim)
                   {
                     shared_ptr<OCCGeometry> geo;
                     if(EndsWith(filename, ".step") || EndsWith(filename, ".stp"))
@@ -145,9 +145,11 @@ DLL_HEADER void ExportNgOCC(py::module &m)
                       geo.reset(LoadOCC_IGES(filename));
                     else
                       throw Exception("Cannot load file " + filename + "\nValid formats are: step, stp, brep, iges");
+                    if(dim<3)
+                      geo->SetDimension(dim);
                     ng_geometry = geo;
                     return geo;
-                  }), py::arg("filename"),
+                  }), py::arg("filename"), py::arg("dim")=3,
         "Load OCC geometry from step, brep or iges file")
     .def(NGSPickle<OCCGeometry>())
     .def("Glue", &OCCGeometry::GlueGeometry)