csg2d interface

libsrc/geom2d/csg2d.cpp

@@ -15,6 +15,22 @@ namespace netgen
 
 constexpr static double EPSILON=0.000000001;
 
+void ComputeWeight( Spline & s, Point<2> p )
+{
+  Point<2> a = s.StartPI();
+  Point<2> b = s.TangentPoint();
+  Point<2> c = s.EndPI();
+
+  double A = (p[1]-a[1])*(b[0]-p[0]) - (p[0]-a[0])*(b[1]-p[1]);
+  double B = (p[1]-c[1])*(b[0]-p[0]) - (p[0]-c[0])*(b[1]-p[1]);
+  double det = sqrt(-A*B);
+  double tt = (B-det)/(A+det);
+  auto v = b-p;
+  int dim = fabs(v[0]) > fabs(v[1]) ? 0 : 1;
+  double weight = fabs(tt*(p[dim]-a[dim])/v[dim] + 1.0/tt*(p[dim]-c[dim])/v[dim]);
+  s.SetWeight(weight);
+}
+
 void ToggleLabel(EntryExitLabel& status)
 {
   if (status == ENTRY)
@@ -61,14 +77,7 @@ Spline Split( const Spline & s, double t0, double t1 )
 
   // compute weight of new spline such that p lies on it
   Point<2> p = s.GetPoint(0.5*(t0+t1));
-  double A = (p[1]-a[1])*(b[0]-p[0]) - (p[0]-a[0])*(b[1]-p[1]);
+  ComputeWeight(res, p);
-  double B = (p[1]-c[1])*(b[0]-p[0]) - (p[0]-c[0])*(b[1]-p[1]);
-  double det = sqrt(-A*B);
-  double tt = (B-det)/(A+det);
-  auto v = b-p;
-  int dim = fabs(v[0]) > fabs(v[1]) ? 0 : 1;
-  double weight = fabs(tt*(p[dim]-a[dim])/v[dim] + 1.0/tt*(p[dim]-c[dim])/v[dim]);
-  res.SetWeight(weight);
   return res;
 }
 
@@ -1222,45 +1231,35 @@ Loop RectanglePoly(double x0, double x1, double y0, double y1, string bc)
   return r;
 }
 
-Solid2d Rectangle(double x0, double x1, double y0, double y1, string name, string bc)
+Solid2d Rectangle(Point<2> p0, Point<2> p1, string name, string bc)
 {
-  Solid2d s;
+  using P = Point<2>;
-  s.name = name;
+  return { {p0, P{p1[0],p0[1]}, p1, P{p0[0],p1[1]}}, name, bc };
-  s.polys.Append(RectanglePoly(x0,x1,y0,y1, bc));
-  s.SetBC(bc);
-  return s;
 }
 
-Solid2d Circle(double x, double y, double r, string name, string bc)
+Solid2d Circle(Point<2> center, double r, string name, string bc)
 {
-  Solid2d s;
+  double x = center[0];
-  s.name = name;
+  double y = center[1];
-  Loop poly;
+  using P = Point<2>;
 
-  Point<2> ps[] =
+  Point<2> p[] =
   {
     {x+r, y+0},
-    {x+r, y+r},
     {x+0, y+r},
-    {x-r, y+r},
     {x-r, y+0},
-    {x-r, y-r},
     {x+0, y-r},
-    {x+r, y-r}
   };
 
-  for (auto i : IntRange(4))
+  EdgeInfo cp[] =
   {
-    int i0 = 2*i;
+    P{x+r, y+r},
-    int i1 = (i0+1)%8;
+    P{x-r, y+r},
-    int i2 = (i0+2)%8;
+    P{x-r, y-r},
-    auto & v0 = poly.Append( ps[i0] );
+    P{x+r, y-r}
-    v0.spline = { ps[i0], ps[i1], ps[i2] };
+  };
-  }
 
-  s.polys.Append(poly);
+  return Solid2d( { p[0], cp[0], p[1], cp[1], p[2], cp[2], p[3], cp[3] }, name, bc );
-  s.SetBC(bc);
-  return s;
 }
 
 Solid2d AddIntersectionPoints ( Solid2d s1, Solid2d s2 )
@@ -1322,7 +1321,7 @@ Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect)
   return res;
 }
 
-Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_)
+Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_, string bc)
     : name(name_)
 {
   Loop l;
@@ -1336,6 +1335,9 @@ Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, strin
 
   for(auto v : l.Vertices(ALL))
     {
+      if(v->info.bc==BC_DEFAULT)
+          v->info.bc = bc;
+
       v->bc = v->info.bc;
       if(v->info.control_point)
         {
@@ -1346,7 +1348,7 @@ Solid2d :: Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, strin
   polys.Append(l);
 }
 
-Solid2d Solid2d :: operator+(Solid2d & other)
+Solid2d Solid2d :: operator+(const Solid2d & other) const
 {
   if(polys.Size()==0)
     return other;
@@ -1356,16 +1358,17 @@ Solid2d Solid2d :: operator+(Solid2d & other)
   return res;
 }
 
-Solid2d Solid2d :: operator*(Solid2d & other)
+Solid2d Solid2d :: operator*(const Solid2d & other) const
 {
   auto res = ClipSolids(*this, other, true);
   res.name = name;
   return res;
 }
 
-Solid2d Solid2d :: operator-(Solid2d other)
+Solid2d Solid2d :: operator-(const Solid2d & other_) const
 {
   // TODO: Check dimensions of solids with bounding box
+  Solid2d other = other_;
   other.Append(RectanglePoly(-1e8, 1e8, -1e8, 1e8, "JUST_FOR_CLIPPING"));
   auto res = ClipSolids(*this, other);
 
@@ -1379,6 +1382,54 @@ Solid2d Solid2d :: operator-(Solid2d other)
   return res;
 }
 
+Solid2d Solid2d :: operator+=(const Solid2d & other)
+{
+  *this = *this + other;
+  return *this;
+}
+
+Solid2d Solid2d :: operator*=(const Solid2d & other)
+{
+  *this = *this * other;
+  return *this;
+}
+
+Solid2d Solid2d :: operator-=(const Solid2d & other)
+{
+  *this = *this - other;
+  return *this;
+}
+
+Solid2d & Solid2d :: Move( Vec<2> v )
+{
+  return Transform( [v](Point<2> p) -> Point<2> { return p+v; } );
+}
+
+Solid2d & Solid2d :: Scale( double sx, double sy )
+{
+  if(sy==0.0)
+      sy=sx;
+  return Transform( [sx,sy](Point<2> p) -> Point<2> { return{p[0]*sx, p[1]*sy}; } );
+}
+
+Solid2d & Solid2d :: RotateRad( double ang, Point<2> center )
+{
+  double sina = sin(ang);
+  double cosa = cos(ang);
+  Vec<2> c = { center[0], center[1] };
+  return Transform( [c, sina, cosa](Point<2> p) -> Point<2>
+      {
+          p -= c;
+          double x = p[0];
+          double y = p[1];
+          p[0] = cosa*x+sina*y;
+          p[1] = -sina*x+cosa*y;
+          p += c;
+          return p;
+      } );
+}
+
+
 bool Solid2d :: IsInside( Point<2> r ) const
 {
   int w = 0;
@@ -1406,7 +1457,6 @@ bool Solid2d :: IsRightInside( const Vertex & p0 )
   return IsInside(q);
 }
 
-
 shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
 {
   static Timer t_intersections("CSG2d - AddIntersections()");
@@ -1519,20 +1569,24 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
         auto li = s.IsLeftInside(p0);
         auto ri = s.IsRightInside(p0);
 
+        auto & ls = seg_map[{pi0,pi1,pi2}];
+        ls.p0 = pi0;
+        ls.p1 = pi1;
+        ls.p2 = pi2;
+        ls.weight = weight;
+
+        if(bcmap.count(p0.bc)==0)
+          bcmap[p0.bc] = bcmap.size()+1;
+
+        if(ls.bc==0 || p0.bc != BC_DEFAULT)
+            ls.bc = bcmap[p0.bc];
+
         if(li!=ri)
         {
-          auto & ls = seg_map[{pi0,pi1,pi2}];
-          ls.p0 = pi0;
-          ls.p1 = pi1;
-          ls.p2 = pi2;
-          ls.weight = weight;
           if(s.IsLeftInside(p0) == flip)
             ls.left = dom;
           else
             ls.right = dom;
-          if(bcmap.count(p0.bc)==0)
-            bcmap[p0.bc] = bcmap.size()+1;
-          ls.bc = bcmap[p0.bc];
         }
       }
     }
@@ -1570,4 +1624,13 @@ shared_ptr<netgen::SplineGeometry2d> CSG2d :: GenerateSplineGeometry()
   }
   return geo;
 }
+
+shared_ptr<netgen::Mesh> CSG2d :: GenerateMesh(MeshingParameters & mp)
+{
+  auto geo = GenerateSplineGeometry();
+  auto mesh = make_shared<netgen::Mesh>();
+  geo->GenerateMesh(mesh, mp);
+  return mesh;
+}
+
 }

libsrc/geom2d/csg2d.hpp

@@ -11,12 +11,15 @@ namespace netgen
 using namespace ngcore;
 using netgen::Point;
 using netgen::Vec;
+using Spline = SplineSeg3<2>;
 
 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]);
 }
 
+// compute weight of spline such that p lies on it
+void ComputeWeight( Spline & s, Point<2> p );
 
 enum IntersectionType
 {                 // types of intersection (detected in the first phase)
@@ -89,8 +92,6 @@ struct EdgeInfo
     }
 };
 
-using Spline = SplineSeg3<2>;
-
 struct Vertex : Point<2>
 {
   Vertex (Point<2> p) : Point<2>(p) {}
@@ -565,15 +566,20 @@ struct Solid2d
 {
   Array<Loop> polys;
 
-  string name = "";
+  string name = MAT_DEFAULT;
 
   Solid2d() = default;
   Solid2d(string name_) : name(name_) {}
-  Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_);
+  Solid2d(const Array<std::variant<Point<2>, EdgeInfo>> & points, string name_=MAT_DEFAULT, string bc_=BC_DEFAULT);
 
-  Solid2d operator+(Solid2d & other);
+  Solid2d operator+(const Solid2d & other) const;
-  Solid2d operator*(Solid2d & other);
+  Solid2d operator*(const Solid2d & other) const;
-  Solid2d operator-(Solid2d other);
+  Solid2d operator-(const Solid2d & other) const;
+
+  Solid2d& operator=(const Solid2d & other) = default;
+  Solid2d operator+=(const Solid2d & other);
+  Solid2d operator*=(const Solid2d & other);
+  Solid2d operator-=(const Solid2d & other);
 
   void Append( const Loop & poly )
   {
@@ -584,14 +590,36 @@ struct Solid2d
   bool IsLeftInside( const Vertex & p0 );
   bool IsRightInside( const Vertex & p0 );
 
-  void SetBC(string bc)
+  template<typename TFunc>
-  {
+  Solid2d & Transform( const TFunc & func )
-    for(auto & p : polys)
+    {
-      for(auto v : p.Vertices(ALL))
+      for(auto & poly : polys)
-        v->bc = bc;
+          for(auto v : poly.Vertices(ALL))
-  }
+            {
+              auto p = func(*v);
+              (*v)[0] = p[0];
+              (*v)[1] = p[1];
+              if(v->spline)
+                {
+                  auto &s = *v->spline;
+                  auto pmid = func(s.GetPoint(0.5));
+                  s = Spline(func(s.StartPI()), func(s.TangentPoint()), func(s.EndPI()));
+                  ComputeWeight(s, pmid);
+                }
+            }
+      return *this;
+    }
+
+  Solid2d & Move( Vec<2> v );
+  Solid2d & Scale( double sx, double sy=0.0 );
+  Solid2d & RotateRad( double ang, Point<2> center = {0,0} );
+  Solid2d & RotateDeg( double ang, Point<2> center = {0,0} )
+    {
+      return RotateRad( ang/180.*M_PI );
+    }
 };
 
+
 class CSG2d
 {
   public:
@@ -603,10 +631,11 @@ class CSG2d
     }
 
     shared_ptr<netgen::SplineGeometry2d> GenerateSplineGeometry();
+    shared_ptr<netgen::Mesh> GenerateMesh(MeshingParameters & mp);
 };
 
-Solid2d Circle(double x, double y, double r, string name="", string bc="");
+Solid2d Circle( Point<2> center, double r, string name="", string bc="");
-Solid2d Rectangle(double x0, double x1, double y0, double y1, string name, string bc);
+Solid2d Rectangle( Point<2> p0, Point<2> p1, string mat=MAT_DEFAULT, string bc=BC_DEFAULT );
 
 Solid2d AddIntersectionPoints ( Solid2d s1, Solid2d s2 );
 Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect=true );

libsrc/geom2d/python_geom2d.cpp

@@ -399,21 +399,30 @@ DLL_HEADER void ExportGeom2d(py::module &m)
   
   py::class_<Solid2d>(m, "Solid2d")
     .def(py::init<>())
-    .def(py::init<Array<std::variant<Point<2>, EdgeInfo>>, std::string>(), py::arg("points"), py::arg("mat")=MAT_DEFAULT)
+    .def(py::init<Array<std::variant<Point<2>, EdgeInfo>>, std::string, std::string>(), py::arg("points"), py::arg("mat")=MAT_DEFAULT, py::arg("bc")=BC_DEFAULT)
     .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("Copy", [](Solid2d & self) -> Solid2d { return self; })
-    .def("SetBC", &Solid2d::SetBC)
+    .def("Move", &Solid2d::Move)
+    .def("Scale", &Solid2d::Scale)
+    .def("Rotate", [](Solid2d & self, optional<double> rad, optional<double> deg, Point<2> center )
+        {
+          if(rad)
+            self.RotateRad(*rad, center);
+          if(deg)
+            self.RotateDeg(*deg, center);
+          return self;
+        }, py::arg("rad")=nullopt, py::arg("deg")=nullopt, py::arg("center")=Point<2>{0,0})
     ;
   
 
-  m.def("Rectangle", [](double x0, double x1, double y0, double y1, string bc, string mat)
+  m.def("Rectangle", [](Point<2> pmin, Point<2> pmax, string bc, string mat)
-		  { return Rectangle(x0,x1,y0,y1,bc,mat); },
+		  { return Rectangle(pmin, pmax, bc,mat); },
-		  py::arg("x0"), py::arg("x1"), py::arg("y0"), py::arg("y1"), py::arg("bc")="", py::arg("mat")=""
+		  py::arg("pmin"), py::arg("pmax"), py::arg("bc")=BC_DEFAULT, py::arg("mat")=MAT_DEFAULT
        );
-  m.def("Circle", Circle, py::arg("x"), py::arg("y"), py::arg("r"), py::arg("bc")="", py::arg("mat")="");
+  m.def("Circle", Circle, py::arg("center"), py::arg("radius"), py::arg("bc")=BC_DEFAULT, py::arg("mat")=MAT_DEFAULT);
 
   py::class_<CSG2d>(m, "CSG2d")
     .def(py::init<>())