Merge branch 'csg2d' into 'master'

CSG for 2D

See merge request jschoeberl/netgen!332

libsrc/geom2d/CMakeLists.txt

@@ -1,5 +1,5 @@
 add_definitions(-DNGLIB_EXPORTS)
-add_library(geom2d ${NG_LIB_TYPE} genmesh2d.cpp geometry2d.cpp python_geom2d.cpp )
+add_library(geom2d ${NG_LIB_TYPE} csg2d.cpp genmesh2d.cpp geometry2d.cpp python_geom2d.cpp )
 if(APPLE)
     set_target_properties( geom2d PROPERTIES SUFFIX ".so")
 endif(APPLE)
@@ -20,6 +20,6 @@ endif(USE_GUI)
 
 install(FILES
   geometry2d.hpp spline2d.hpp
-  vsgeom2d.hpp
+  vsgeom2d.hpp csg2d.hpp
   DESTINATION ${NG_INSTALL_DIR_INCLUDE}/geom2d COMPONENT netgen_devel
 )

libsrc/geom2d/csg2d.hpp

@@ -0,0 +1,583 @@
+#ifndef NETGEN_CSG2D_HPP_INCLUDED
+#define NETGEN_CSG2D_HPP_INCLUDED
+
+#include "geometry2d.hpp"
+
+namespace netgen
+{
+
+using namespace ngcore;
+using netgen::Point;
+using netgen::Vec;
+
+inline double Area(const Point<2>& P, const Point<2>& Q, const Point<2>& R)
+{
+  return (Q[0]-P[0]) * (R[1]-P[1]) - (Q[1]-P[1]) * (R[0]-P[0]);
+}
+
+
+enum IntersectionType
+{                 // types of intersection (detected in the first phase)
+  NO_INTERSECTION = 0,
+  X_INTERSECTION,
+  T_INTERSECTION_Q,
+  T_INTERSECTION_P,
+  V_INTERSECTION,
+  X_OVERLAP,
+  T_OVERLAP_Q,
+  T_OVERLAP_P,
+  V_OVERLAP
+};
+
+enum IntersectionLabel
+{      // for the classification of intersection vertices in the second phase
+  NONE,
+  CROSSING,
+  BOUNCING,
+  LEFT_ON,
+  RIGHT_ON,
+  ON_ON,
+  ON_LEFT,
+  ON_RIGHT,
+  DELAYED_CROSSING,
+  DELAYED_BOUNCING
+};
+
+enum EntryExitLabel
+{         // for marking intersection vertices as "entry" or "exit"
+  EXIT,
+  ENTRY,
+  NEITHER
+};
+
+enum IteratorType
+{
+  SOURCE,
+  INTERSECTION,
+  CROSSING_INTERSECTION,
+  ALL
+};
+
+
+using Spline = SplineSeg3<2>;
+
+struct Vertex : Point<2>
+{
+  Vertex (Point<2> p) : Point<2>(p) {}
+
+  Vertex * prev = nullptr;
+  Vertex * next = nullptr;
+  unique_ptr<Vertex> pnext = nullptr;
+  Vertex * neighbour = nullptr; // same vertex in other polygon (at intersections)
+  double lam = -1.0;
+  bool is_intersection = false;
+  bool is_source = false;
+
+  IntersectionLabel label = NONE;    // type of intersection vertex
+  EntryExitLabel enex = NEITHER;        // entry/exit "flag"
+
+  string bc = "";
+
+  // In case the edge this - next is curved, store the spline information here
+  optional<Spline> spline = nullopt;
+
+  Vertex * Insert(Point<2> p, double lam = -1.0);
+
+  void Link( Vertex * v )
+  {
+    neighbour = v;
+    v->neighbour = this;
+    is_intersection = true;
+    v->is_intersection = true;
+  }
+};
+
+struct VertexIterator
+{
+  struct iterator
+  {
+    iterator(Vertex* root, IteratorType IterType) :
+      root(root), V(NULL), iterType(IterType)
+    {
+      if (root == NULL)
+        return;
+
+      if (nextVertex() == NULL)         // no (source/intersection) vertex found
+        root = V = NULL;                // -> mark iterator as "end"
+    }
+
+    const iterator& operator++()
+    {
+      nextVertex();
+      return *this;
+    }
+
+    Vertex* operator*()
+    {
+      return V;
+    }
+
+    bool operator!=(const iterator& other) const
+    {
+      return (root != other.root) || (V != other.V);
+    }
+
+    private:
+    Vertex* root;
+    Vertex* V;
+    IteratorType iterType;
+
+    //
+    // find the next vertex
+    // if iterType is ALL, then it is just the next vertex
+    // if iterType is SOURCE, then it is the next source vertex
+    // if iterType is INTERSECTION, then it is the next intersection vertex
+    // if iterType is CROSSING_INTERSECTION, then it is the next intersection vertex with CROSSING label
+    //
+    Vertex* nextVertex()
+    {
+      bool nextFound = false;
+
+      if (V == NULL)
+      {                  // find first (source/intersection) vertex
+        V = root;
+        switch(iterType)
+        {
+          case ALL:
+            nextFound = true;
+            break;
+          case SOURCE:
+            if (V->is_source)
+              nextFound = true;
+            break;
+          case INTERSECTION:
+            if (V->is_intersection)
+              nextFound = true;
+            break;
+          case CROSSING_INTERSECTION:
+            if (V->is_intersection && (V->label == CROSSING))
+              nextFound = true;
+            break;
+        }
+      }
+
+      while (!nextFound)
+      {              // find next (source/intersection) vertex
+        switch(iterType)
+        {
+          case ALL:
+            V = V->next;
+            break;
+          case SOURCE:
+            do {
+              V = V->next;
+            } while (!V->is_source && V != root);
+            break;
+          case INTERSECTION:
+            do {
+              V = V->next;
+            } while (!V->is_intersection && V != root);
+            break;
+          case CROSSING_INTERSECTION:
+            do {
+              V = V->next;
+            } while ( ( !V->is_intersection || (V->label != CROSSING) ) && V != root);
+            break;
+        }
+
+        if (V == root)
+        {                // back at the root vertex?
+          root = V = NULL;              // -> mark iterator as "end"
+          return(V);
+        }
+
+        switch(iterType)
+        {
+          case ALL:
+            nextFound = true;
+            break;
+          case SOURCE:
+            if (V->is_source)
+              nextFound = true;
+            break;
+          case INTERSECTION:
+            if (V->is_intersection)
+              nextFound = true;
+            break;
+          case CROSSING_INTERSECTION:
+            if (V->is_intersection && (V->label == CROSSING))
+              nextFound = true;
+            break;
+        }
+      }
+      return(V);
+    }
+  };
+
+  public:
+  VertexIterator() : root(NULL) {};
+
+  iterator begin() { return iterator(root, iterType); }
+  iterator end()   { return iterator(NULL, iterType); }
+
+  Vertex* root;
+  IteratorType iterType;
+};
+
+
+struct Edge
+{
+  Vertex * v0 = nullptr;
+  Vertex * v1 = nullptr;
+
+  Edge (Vertex* v, Vertex* w) : v0(v), v1(w) { };
+};
+
+struct EdgeIterator
+{
+  struct iterator
+  {
+    iterator(Vertex* root, IteratorType IterType) :
+      root(root), one(NULL), two(NULL), iterType(IterType)
+    {
+      if (root == NULL)
+        return;
+
+      if (nextEdge() == NULL)           // no source edge found
+        root = one = two = NULL;        // -> mark iterator as "end"
+    }
+
+    const iterator& operator++() { nextEdge(); return *this; }
+
+    Edge operator*()
+    {
+      return Edge(one,two);
+    }
+
+    bool operator!=(const iterator& other) const
+    {
+      return (root != other.root) || (one != other.one) || (two != other.two);
+    }
+
+    private:
+    Vertex* root;
+    Vertex* one;
+    Vertex* two;
+    IteratorType iterType;
+
+    //
+    // find the next vertex, starting at curr
+    // if iterType is ALL, then it is just the next vertex
+    // if iterType is SOURCE, then it is the next source vertex
+    //
+    Vertex* nextVertex(Vertex* curr)
+    {
+      if (curr == NULL)
+        return(NULL);
+
+      switch(iterType)
+      {
+        case ALL:
+          curr = curr->next;
+          break;
+
+        case SOURCE:
+          do {
+            curr = curr->next;
+          } while (!curr->is_source);
+          break;
+        default:
+          ;
+      }
+
+      return(curr);
+    }
+
+    //
+    // find the next edge
+    //
+    Vertex* nextEdge()
+    {
+      if (root == NULL)                 // empty polygon?
+        return (NULL);
+
+      if (one == NULL)
+      {                // find one (source) vertex
+        one = root;                     // note: root is always a (source) vertex
+        two = nextVertex(one);
+        if (two == one)                 // just one (source) vertex
+          return(NULL);                 // -> no (source) edges
+        return(one);
+      }
+
+      if (two == root)
+      {                // back at the root vertex?
+        root = one = two = NULL;        // -> mark iterator as "end"
+        return(NULL);
+      }
+
+      one = two;
+      two = nextVertex(one);
+
+      return (one);
+    }
+  };
+
+  public:
+  EdgeIterator() : root(NULL) {};
+
+  iterator begin() { return iterator(root, iterType); }
+  iterator end()   { return iterator(NULL, iterType); }
+
+  Vertex* root;
+  IteratorType iterType;
+};
+
+
+inline int CalcSide( const Point<2> & p0, const Point<2> & p1, const Point<2> & r )
+{
+  if ( (p0[1] < r[1]) != (p1[1] < r[1]) )
+  {
+    if (p0[0] >= r[0])
+    {
+      if (p1[0] > r[0])
+        return 2 * (p1[1] > p0[1]) - 1;
+      else
+        if ( (Area(p0,p1,r) > 0) == (p1[1] > p0[1]) )
+          return 2 * (p1[1] > p0[1]) - 1;
+    }
+    else
+    {
+      if (p1[0] > r[0])
+        if ( (Area(p0,p1,r) > 0) == (p1[1] > p0[1]) )
+          return 2 * (p1[1] > p0[1]) - 1;
+    }
+  }
+  return 0;
+}
+
+struct Polygon2d
+{
+  unique_ptr<Vertex> first = nullptr;
+
+  Polygon2d() = default;
+
+  Polygon2d(const Polygon2d & p)
+    : first(nullptr)
+  {
+    for(auto v : p.Vertices(ALL))
+      AppendVertex(*v);
+  }
+
+  Polygon2d & operator=(const Polygon2d & p)
+  {
+    first = nullptr;
+    if(p.first)
+      for(const auto v : p.Vertices(ALL))
+        AppendVertex(*v);
+    return *this;
+  }
+
+  Vertex & AppendVertex(const Vertex & v)
+  {
+    auto & vnew = Append( static_cast<Point<2>>(v), true );
+    vnew.bc = v.bc;
+    if(v.spline)
+      vnew.spline = *v.spline;
+    return vnew;
+  }
+
+  Vertex & Append(Point<2> p, bool source = false)
+  {
+    Vertex * vnew;
+    if(first==nullptr)
+    {
+      first = make_unique<Vertex>(p);
+      first->next = first.get();
+      first->prev = first.get();
+      vnew = first.get();
+    }
+    else
+    {
+      vnew = first->prev->Insert(p);
+    }
+
+    vnew->is_source = source;
+    //     cout << "size after " << Size() << endl;
+    return *vnew;
+  }
+
+  void Remove (Vertex* v)
+  {
+    v->prev->next = v->next;
+    v->next->prev = v->prev;
+    if(first.get() == v)
+      first = std::move(v->pnext);
+    else
+      v->prev->pnext = std::move(v->pnext);
+  }
+
+  bool IsInside( Point<2> r ) const
+  {
+    int w = 0;
+    for(auto e : Edges(ALL))
+      w += CalcSide(*e.v0, *e.v1, r);
+    return ( (w % 2) != 0 );
+  }
+
+  EdgeIterator Edges(IteratorType iterType) const
+  {
+    EdgeIterator it;
+    it.iterType = iterType;
+    it.root = first.get();
+    return it;
+  }
+
+  VertexIterator Vertices(IteratorType iterType, Vertex* first_ = nullptr) const
+  {
+    VertexIterator it;
+    it.iterType = iterType;
+    it.root = (first_ == nullptr) ? first.get() : first_;
+    return it;
+  }
+
+  //
+  // check, if all vertices have the ON_ON label
+  //
+  bool allOnOn()
+  {
+    for (Vertex* v : Vertices(ALL))
+      if (v->label != ON_ON)
+        return(false);
+    return(true);
+  }
+
+  //
+  // check, if the polygon does not contain any crossing intersection vertex
+  // or crossing intersection chain or (if we want to compute the union instead
+  // of the intersection) a bouncing vertex or a bouncing intersection chain
+  //
+  bool noCrossingVertex(bool union_case = false)
+  {
+    for (Vertex* v : Vertices(ALL))
+      if (v->is_intersection)
+      {
+        if ( (v->label == CROSSING) || (v->label == DELAYED_CROSSING) )
+          return(false);
+
+        if (union_case && ( (v->label == BOUNCING) || (v->label == DELAYED_BOUNCING) ) )
+          return(false);
+      }
+    return(true);
+  }
+
+  //
+  // return a non-intersection point
+  //
+  Point<2> getNonIntersectionPoint()
+  {
+    for (Vertex* v : Vertices(ALL))
+      if (!v->is_intersection)
+        return *v;
+
+    // no non-intersection vertex found -> find suitable edge midpoint
+    for (Vertex* v : Vertices(ALL))
+      // make sure that edge from V to V->next is not collinear with other polygon
+      if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) )
+        // return edge midpoint
+        return Center(*v, *v->next);
+    throw Exception("no point found");
+  }
+
+  //
+  // return and insert a non-intersection vertex
+  //
+  Vertex* getNonIntersectionVertex()
+  {
+    for (Vertex* v : Vertices(ALL))
+      if (!v->is_intersection)
+        return(v);
+
+    // no non-intersection vertex found -> generate and return temporary vertex
+    for (Vertex* v : Vertices(ALL))
+      // make sure that edge from V to V->next is not collinear with other polygon
+      if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) )
+      {
+        // add edge midpoint as temporary vertex
+        auto p = Center(*v, *v->next);
+        return v->Insert(p);
+      }
+    return(NULL);
+  }
+
+  void SetBC(string bc)
+  {
+    for(auto v : Vertices(ALL))
+      v->bc = bc;
+  }
+
+  size_t Size() const
+  {
+    if(first==nullptr) return 0;
+
+    size_t cnt = 0;
+
+    for(auto v : Vertices(ALL))
+      cnt++;
+
+    return cnt;
+  }
+};
+
+
+struct Solid2d
+{
+  Array<Polygon2d> polys;
+
+  string name = "";
+
+  Solid2d() = default;
+  Solid2d(string name_) : name(name_) {}
+
+  Solid2d operator+(Solid2d & other);
+  Solid2d operator*(Solid2d & other);
+  Solid2d operator-(Solid2d other);
+
+  void Append( const Polygon2d & poly )
+  {
+    polys.Append(poly);
+  }
+
+  bool IsInside( Point<2> r ) const;
+  bool IsLeftInside( const Vertex & p0 );
+  bool IsRightInside( const Vertex & p0 );
+
+  void SetBC(string bc)
+  {
+    for(auto & p : polys)
+      for(auto v : p.Vertices(ALL))
+        v->bc = bc;
+  }
+};
+
+class CSG2d
+{
+  public:
+    Array<Solid2d> solids;
+
+    void Add ( Solid2d s )
+    {
+      solids.Append(s);
+    }
+
+    shared_ptr<netgen::SplineGeometry2d> GenerateSplineGeometry();
+};
+
+Solid2d Circle(double x, double y, double r, string name="", string bc="");
+Solid2d Rectangle(double x0, double x1, double y0, double y1, string name, string bc);
+
+Solid2d AddIntersectionPoints ( Solid2d s1, Solid2d s2 );
+Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect=true );
+
+}
+#endif // NETGEN_CSG2D_HPP_INCLUDED

libsrc/geom2d/geometry2d.hpp

@@ -9,6 +9,7 @@
 
 #include <myadt.hpp>
 #include <gprim.hpp>
+#include <meshing.hpp>
 
 
 // #include "../gprim/spline.hpp"

libsrc/geom2d/python_geom2d.cpp

@@ -6,6 +6,7 @@
 
 #include <meshing.hpp>
 #include <geometry2d.hpp>
+#include <csg2d.hpp>
 
 using namespace netgen;
 
@@ -265,6 +266,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
                  {
                    double len = self.splines[i]->Length();
                    int n = floor(len/(0.05*min(xdist,ydist)));
+                   n = max(3, n);
                    lst.push_back(self.splines[i]->StartPI());
                    for (int j = 1; j < n; j++){
                      lst.push_back(self.splines[i]->GetPoint(j*1./n));
@@ -395,6 +397,39 @@ DLL_HEADER void ExportGeom2d(py::module &m)
     .def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing)
     ;
   
+  py::class_<Solid2d>(m, "Solid2d")
+    .def(py::init<>())
+    .def(py::init<std::string>())
+    .def_readwrite("name", &Solid2d::name)
+    .def("__mul__", [](Solid2d & self, Solid2d & other) { return self*other; })
+    .def("__add__", [](Solid2d & self, Solid2d & other) { return self+other; })
+    .def("__sub__", [](Solid2d & self, Solid2d & other) { return self-other; })
+    .def("Append", &Solid2d::Append)
+    .def("SetBC", &Solid2d::SetBC)
+    ;
+  
+  py::class_<Polygon2d>(m, "Polygon2d")
+    .def(py::init<>())
+    .def("SetBC", &Polygon2d::SetBC)
+    .def("Append", [](Polygon2d & self, double x, double y)
+        {
+        self.Append({x,y});
+        })
+    ;
+
+
+  m.def("Rectangle", [](double x0, double x1, double y0, double y1, string bc, string mat)
+		  { return Rectangle(x0,x1,y0,y1,bc,mat); },
+		  py::arg("x0"), py::arg("x1"), py::arg("y0"), py::arg("y1"), py::arg("bc")="", py::arg("mat")=""
+       );
+  m.def("Circle", Circle, py::arg("x"), py::arg("y"), py::arg("r"), py::arg("bc")="", py::arg("mat")="");
+
+  py::class_<CSG2d>(m, "CSG2d")
+    .def(py::init<>())
+    .def("GenerateSplineGeometry", &CSG2d::GenerateSplineGeometry)
+    .def("Add", &CSG2d::Add)
+    ;
+
 }
 
 PYBIND11_MODULE(libgeom2d, m) {

libsrc/gprim/geom2d.cpp

@@ -251,239 +251,4 @@ int PTRIANGLE2D :: IsIn (const Point2d & p) const
 }
 
 #endif
-
-
-
-
-
-
-Polygon2d :: Polygon2d ()
-{
-  ;
-}
-
-Polygon2d :: ~Polygon2d ()
-{
-  ;
-}
-
-void Polygon2d :: AddPoint (const Point2d & p)
-{ 
-  points.Append(p); 
-}
-
-
-double Polygon2d :: HArea () const
-{
-  int i;
-  double ar = 0;
-  for (i = 1; i <= points.Size(); i++)
-    {
-      const Point2d & p1 = points.Get(i);
-      const Point2d & p2 = points.Get(i%points.Size()+1);
-      ar += 
-	(p2.X()-p1.X()) * p1.Y() -
-	(p2.Y()-p1.Y()) * p1.X();
-    }
-  return ar/2;
-  /*
-  CURSOR c;
-  double ar = 0;
-  Point2d * p1, * p2, p0 = Point2d(0, 0);
-  Vec2d v1, v2 = Vec2d(1, 0);
-
-  p2 = points[points.Last()];
-  for (c = points.First(); c != points.Head(); c++)
-    {
-    p1 = p2;
-    p2 = points[c];
-    ar += Cross ( (*p2-*p1), (*p1 - p0));
-    }
-  return ar / 2;
-  */
-}
-
-
-int Polygon2d :: IsOn (const Point2d & p) const
-{
-  int i;
-  for (i = 1; i <= points.Size(); i++)
-    {
-      const Point2d & p1 = points.Get(i);
-      const Point2d & p2 = points.Get(i%points.Size()+1);
-      if (IsOnLine (Line2d(p1, p2), p)) return 1;
-    }
-  return 0;
-  /*
-  CURSOR c;
-  Point2d * p1, * p2;
-  
-  p2 = points[points.Last()];
-  for (c = points.First(); c != points.Head(); c++)
-    {
-      p1 = p2;
-      p2 = points[c];
-      if (IsOnLine (Line2d(*p1, *p2), p)) return 1;
-    }
-  return 0;
-  */
-}
-
-
-int Polygon2d :: IsIn (const Point2d & p) const
-{
-  int i;
-  double sum = 0, ang;
-  for (i = 1; i <= points.Size(); i++)
-    {
-      const Point2d & p1 = points.Get(i);
-      const Point2d & p2 = points.Get(i%points.Size()+1);
-      ang = Angle ( (p1 - p), (p2 - p) );
-      if (ang > M_PI) ang -= 2 * M_PI;
-      sum += ang;
-    }
-  return fabs(sum) > M_PI;
-  /*
-  CURSOR c;
-  Point2d * p1, * p2;
-  double sum = 0, ang;
-
-  p2 = points[points.Last()];
-  for (c = points.First(); c != points.Head(); c++)
-    {
-    p1 = p2;
-    p2 = points[c];
-    ang = Angle ( (*p1 - p), (*p2 - p) );
-    if (ang > M_PI) ang -= 2 * M_PI;
-    sum += ang;
-    }
-
-  return fabs(sum) > M_PI;
-  */
-}
-
-int Polygon2d :: IsConvex () const
-  {
-    /*
-  Point2d *p, *pold, *pnew;
-  char cw;
-  CURSOR c;
-
-  if (points.Length() < 3) return 0;
-
-  c = points.Last();
-  p = points[c];
-  c--;
-  pold = points[c];
-  pnew = points[points.First()];
-  cw = ::CW (*pold, *p, *pnew);
-
-  for (c = points.First(); c != points.Head(); c++)
-    {
-    pnew = points[c];
-    if (cw != ::CW (*pold, *p, *pnew))
-      return 0;
-    pold = p;
-    p = pnew;
-    }
-    */
-    return 0;
-  }
-
-
-int Polygon2d :: IsStarPoint (const Point2d & p) const
-  {
-    /*
-  Point2d *pnew, *pold;
-  char cw;
-  CURSOR c;
-
-  if (points.Length() < 3) return 0;
-
-  pold = points[points.Last()];
-  pnew = points[points.First()];
-
-  cw = ::CW (p, *pold, *pnew);
-
-  for (c = points.First(); c != points.Head(); c++)
-    {
-    pnew = points[c];
-    if (cw != ::CW (p, *pold, *pnew))
-      return 0;
-    pold = pnew;
-    }
-  return 1;
-    */
-    return 0;
-  }
-
-
-Point2d Polygon2d :: Center () const
-  {
-    /*
-  double ai, a = 0, x = 0, y = 0;
-  Point2d * p, *p2;
-  Point2d p0 = Point2d(0, 0);
-  CURSOR c;
-
-  p2 = points[points.Last()];
-
-  for (c = points.First(); c != points.Head(); c++)
-    {
-    p = points[c];
-    ai = Cross (*p2 - p0, *p - p0);
-    x += ai / 3 * (p2->X() + p->X());
-    y += ai / 3 * (p2->Y() + p->Y());
-    a+= ai;
-    p2 = p;
-    }
-  if (a != 0)
-    return Point2d (x / a, y / a);
-  else
-    return Point2d (0, 0);
-    */
-    return Point2d (0, 0);
-  }
-
-
-
-Point2d Polygon2d :: EqualAreaPoint () const
-  {
-    /*
-  double a11 = 0, a12 = 0, a21= 0, a22 = 0;
-  double b1 = 0, b2 = 0, dx, dy;
-  double det;
-  Point2d * p, *p2;
-  CURSOR c;
-
-  p = points[points.Last()];
-
-  for (c = points.First(); c != points.Head(); c++)
-    {
-    p2 = p;
-    p = points[c];
-
-    dx = p->X() - p2->X();
-    dy = p->Y() - p2->Y();
-
-    a11 += sqr (dy);
-    a12 -= dx * dy;
-    a21 -= dx * dy;
-    a22 += sqr (dx);
-    b1 -= dy * (p->X() * p2->Y() - p2->X() * p->Y());
-    b2 -= dx * (p->Y() * p2->X() - p2->Y() * p->X());
-    }
-
-  det = a11 * a22 - a21 * a12;
-
-  if (det != 0)
-    return Point2d ( (b1 * a22 - b2 * a12) / det,
-                     (a11 * b2 - a21 * b1) / det);
-  else
-    return Point2d (0, 0);
-*/
-    return Point2d (0, 0);
-  }
-
-
 }

libsrc/gprim/geom2d.hpp

@@ -609,40 +609,6 @@ namespace netgen
 #endif
 
 
-
-  class Polygon2d
-  {
-  protected:
-    NgArray<Point2d> points;
-  
-  public:
-    Polygon2d ();
-    ~Polygon2d ();
-  
-    void AddPoint (const Point2d & p);
-    int GetNP() const { return points.Size(); }
-    void GetPoint (int i, Point2d & p) const
-    { p = points.Get(i); }
-    void GetLine (int i, Point2d & p1, Point2d & p2) const
-    { p1 = points.Get(i); p2 = points.Get(i%points.Size()+1); }
-
-    double Area () const { return fabs (HArea()); }
-    int CW () const { return HArea() > 0; }
-    int CCW () const { return HArea() < 0; }
-  
-    int IsOn (const Point2d & p) const;
-    int IsIn (const Point2d & p) const;
-
-    int IsConvex () const;
-
-    int IsStarPoint (const Point2d & p) const;
-    Point2d Center() const;
-    Point2d EqualAreaPoint () const;
-  private:
-    double HArea () const;
-  };
-
-
   /** Cheap approximation to atan2.
       A monotone function of atan2(x,y) is computed.
   */

libsrc/gprim/spline.hpp

@@ -196,6 +196,10 @@ namespace netgen
     {
       ar & p1 & p2 & p3 & weight & proj_latest_t;
     }
+    ///
+    double GetWeight () const { return weight; }
+    void SetWeight (double w) { weight = w; }
+    ///
     virtual Point<D> GetPoint (double t) const;
     ///
     virtual Vec<D> GetTangent (const double t) const;

py_tutorials/csg2d.py

@@ -0,0 +1,50 @@
+from ngsolve import *
+
+from random import random, seed
+ngsglobals.msg_level = 0
+
+import netgen
+from pyngcore import *
+from netgen.geom2d import *
+
+seed(4)
+
+def GenerateMesh():
+
+    g = CSG2d()
+    g1 = CSG2d()
+    outer = Rectangle(0, 1, 0, 1,"outer","outer")
+    inner = Solid2d()
+
+    for i in range(30):
+        cx = random()
+        cy = random()
+        r  = 0.03+0.05*random()
+        print("Add Circle", i, cx, cy, r, flush = True)
+        circle = Circle(cx, cy, r, "circle"+str(i), "circle"+str(i))
+        g1.Add(circle)
+        inner += circle
+        outer -= circle
+
+
+    g.Add(inner)
+    g.Add(outer)
+    geo = g.GenerateSplineGeometry()
+    Draw(geo)
+
+    # draw this geometry for checking ff the final mesh/geometry is correct
+    # g1.Add(outer)
+    # geo1 = g1.GenerateSplineGeometry()
+    # Draw(geo1)
+
+    print('generate mesh')
+    m = geo.GenerateMesh(maxh=0.1)
+    mesh = Mesh(m)
+    mesh.Curve(3)
+    Draw(mesh)
+
+    return mesh
+
+from ngsolve import Draw
+with PajeTrace():
+    mesh = GenerateMesh()

python/geom2d.py

@@ -1,4 +1,4 @@
-from .libngpy._geom2d import SplineGeometry
+from .libngpy._geom2d import SplineGeometry, Solid2d, Polygon2d, CSG2d, Rectangle, Circle
 from .meshing import meshsize
 
 unit_square = SplineGeometry()