diff --git a/libsrc/geom2d/csg2d.cpp b/libsrc/geom2d/csg2d.cpp index 1f323193..9a1a208c 100644 --- a/libsrc/geom2d/csg2d.cpp +++ b/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]); - 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); + ComputeWeight(res, p); 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; - s.name = name; - s.polys.Append(RectanglePoly(x0,x1,y0,y1, bc)); - s.SetBC(bc); - return s; + using P = Point<2>; + return { {p0, P{p1[0],p0[1]}, p1, P{p0[0],p1[1]}}, name, bc }; } -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; - s.name = name; - Loop poly; + double x = center[0]; + double y = center[1]; + 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; - int i1 = (i0+1)%8; - int i2 = (i0+2)%8; - auto & v0 = poly.Append( ps[i0] ); - v0.spline = { ps[i0], ps[i1], ps[i2] }; - } + P{x+r, y+r}, + P{x-r, y+r}, + P{x-r, y-r}, + P{x+r, y-r} + }; - s.polys.Append(poly); - s.SetBC(bc); - return s; + return Solid2d( { p[0], cp[0], p[1], cp[1], p[2], cp[2], p[3], cp[3] }, name, bc ); } Solid2d AddIntersectionPoints ( Solid2d s1, Solid2d s2 ) @@ -1322,7 +1321,7 @@ Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect) return res; } -Solid2d :: Solid2d(const Array, EdgeInfo>> & points, string name_) +Solid2d :: Solid2d(const Array, EdgeInfo>> & points, string name_, string bc) : name(name_) { Loop l; @@ -1336,6 +1335,9 @@ Solid2d :: Solid2d(const Array, 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, 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 CSG2d :: GenerateSplineGeometry() { static Timer t_intersections("CSG2d - AddIntersections()"); @@ -1519,20 +1569,24 @@ shared_ptr 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 CSG2d :: GenerateSplineGeometry() } return geo; } + +shared_ptr CSG2d :: GenerateMesh(MeshingParameters & mp) +{ + auto geo = GenerateSplineGeometry(); + auto mesh = make_shared(); + geo->GenerateMesh(mesh, mp); + return mesh; +} + } diff --git a/libsrc/geom2d/csg2d.hpp b/libsrc/geom2d/csg2d.hpp index 1994bb27..c56859f8 100644 --- a/libsrc/geom2d/csg2d.hpp +++ b/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 polys; - string name = ""; + string name = MAT_DEFAULT; Solid2d() = default; Solid2d(string name_) : name(name_) {} - Solid2d(const Array, EdgeInfo>> & points, string name_); + Solid2d(const Array, EdgeInfo>> & points, string name_=MAT_DEFAULT, string bc_=BC_DEFAULT); - Solid2d operator+(Solid2d & other); - Solid2d operator*(Solid2d & other); - Solid2d operator-(Solid2d other); + Solid2d operator+(const Solid2d & other) const; + Solid2d operator*(const Solid2d & other) const; + 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) - { - for(auto & p : polys) - for(auto v : p.Vertices(ALL)) - v->bc = bc; - } + template + Solid2d & Transform( const TFunc & func ) + { + for(auto & poly : polys) + 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 GenerateSplineGeometry(); + shared_ptr GenerateMesh(MeshingParameters & mp); }; -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 Circle( Point<2> center, double r, 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 ); diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index fb550d79..59de6c51 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -399,21 +399,30 @@ DLL_HEADER void ExportGeom2d(py::module &m) py::class_(m, "Solid2d") .def(py::init<>()) - .def(py::init, EdgeInfo>>, std::string>(), py::arg("points"), py::arg("mat")=MAT_DEFAULT) + .def(py::init, 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("SetBC", &Solid2d::SetBC) + .def("Copy", [](Solid2d & self) -> Solid2d { return self; }) + .def("Move", &Solid2d::Move) + .def("Scale", &Solid2d::Scale) + .def("Rotate", [](Solid2d & self, optional rad, optional 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) - { 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("Rectangle", [](Point<2> pmin, Point<2> pmax, string bc, string mat) + { return Rectangle(pmin, pmax, bc,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_(m, "CSG2d") .def(py::init<>())