From f2c6a0f8c0a2c7f7946bc4e3b0ec43a8e9b37518 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Mon, 22 Nov 2021 15:23:34 +0100 Subject: [PATCH 01/20] extended wrappers for spline approximation; enable spline surface interpolation --- libsrc/occ/python_occ_shapes.cpp | 252 +++++++++++++++++++++++++++++-- 1 file changed, 243 insertions(+), 9 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index f99bd6e0..f610fdbf 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -35,6 +35,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -685,8 +689,6 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) .value("SHAPE", TopAbs_SHAPE) .export_values() ; - - py::class_ (m, "TopoDS_Shape") @@ -1712,6 +1714,22 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) }) ; + + py::enum_(m, "ShapeContinuity", "Wrapper for OCC enum GeomAbs_Shape") + .value("C0", GeomAbs_Shape::GeomAbs_C0) + .value("C1", GeomAbs_Shape::GeomAbs_C1) + .value("C2", GeomAbs_Shape::GeomAbs_C2) + .value("C3", GeomAbs_Shape::GeomAbs_C3) + .value("CN", GeomAbs_Shape::GeomAbs_CN) + .value("G1", GeomAbs_Shape::GeomAbs_G1) + .value("G2", GeomAbs_Shape::GeomAbs_G2); + + py::enum_(m, "ApproxParamType", "Wrapper for Approx_ParametrizationType") + .value("Centripetal", Approx_ParametrizationType::Approx_Centripetal) + .value("ChordLength", Approx_ParametrizationType::Approx_ChordLength) + .value("IsoParametric", Approx_ParametrizationType::Approx_IsoParametric); + + m.def("HalfSpace", [] (gp_Pnt p, gp_Vec n) { gp_Pln plane(p, n); @@ -1994,14 +2012,230 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return BRepBuilderAPI_MakeEdge(curve).Edge(); }, py::arg("points"), "create Bezier curve"); - m.def("SplineApproximation", [](std::vector pnts, double tol) { - TColgp_Array1OfPnt points(0, pnts.size()-1); - for (int i = 0; i < pnts.size(); i++) - points.SetValue(i, pnts[i]); - GeomAPI_PointsToBSpline builder(points); + m.def("SplineApproximation", [](py::object pnts, Approx_ParametrizationType approx_type, int deg_min, + int deg_max, GeomAbs_Shape continuity, double tol) { + TColgp_Array1OfPnt points(0, 0); + if (py::extract>(pnts).check()) { + std::vector pnt_list{py::extract>(pnts)()}; + points.Resize(0, pnt_list.size()-1, true); + for (int i = 0; i < pnt_list.size(); i++) + points.SetValue(i, pnt_list[i]); + } else if (py::extract>(pnts).check()) { + py::array_t pnt_array{py::extract>(pnts)()}; + if (pnt_array.ndim() != 2) + throw Exception("`points` array must have dimension 2."); + if (pnt_array.shape(1) != 3) + throw Exception("The second dimension must have size 3."); + + points.Resize(0, pnt_array.shape(0)-1, true); + auto pnts_unchecked = pnt_array.unchecked<2>(); + for (int i = 0; i < pnt_array.shape(0); ++i) + points.SetValue(i, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); + } else + throw Exception("Not able to process the data type of points"); + + GeomAPI_PointsToBSpline builder(points, approx_type, deg_min, deg_max, continuity, tol); return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); - }, py::arg("points"), py::arg("tol"), - "Generate spline-curve approximating list of points up to tolerance tol"); + }, + py::arg("points"), + py::arg("approx_type") = Approx_ParametrizationType::Approx_ChordLength, + py::arg("deg_min") = 3, + py::arg("deg_max") = 8, + py::arg("continuity") = GeomAbs_Shape::GeomAbs_C2, + py::arg("tol")=1e-8, + R"delimiter( +Generate a piecewise continuous spline-curve approximating a list of points. + +Parameters +---------- + +points : List[gp_Pnt] or Tuple[gp_Pnt] or np.ndarray[double] + List (or tuple) of gp_Pnt. If a numpy array is provided instead, the data must contain the coordinates + +approx_type : ApproxParamType + Assumption on location of parameters wrt points. + +deg_min : int + Minimum polynomial degree of splines + +deg_max : int + Maxmium polynomial degree of splines + +continuity : ShapeContinuity + Continuity requirement on the approximating surface + +tol : float + Tolerance for the distance from individual points to the approximating curve. + +)delimiter"); + +// CRASHES BADLY during call of Curve()! +// m.def("SplineInterpolation", [](py::object pnts, bool periodic, double tol) { +// cout << "enter" << endl; +// TColgp_Array1OfPnt points(0, 0); +// cout << "points array exists" << endl; +// if (py::extract>(pnts).check()) { +// std::vector pnt_list{py::extract>(pnts)()}; +// points.Resize(1, pnt_list.size(), true); +// for (int i = 0; i < pnt_list.size(); i++) +// points.SetValue(i+1, pnt_list[i]); +// } else if (py::extract>(pnts).check()) { +// py::array_t pnt_array{py::extract>(pnts)()}; +// if (pnt_array.ndim() != 2) +// throw Exception("`points` array must have dimension 2."); +// if (pnt_array.shape(1) != 3) +// throw Exception("The second dimension must have size 3."); +// cout << "resize" << endl; +// points.Resize(1, pnt_array.shape(0), true); +// auto pnts_unchecked = pnt_array.unchecked<2>(); +// for (int i = 0; i < pnt_array.shape(0); ++i) +// points.SetValue(i+1, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); +// cout << "values set" << endl; +// } else +// throw Exception("Not able to process the data type of points"); +// +// TColgp_HArray1OfPnt hpoints{points}; +// const auto _handle = opencascade::handle{&hpoints}; +// cout << "build" << endl; +// GeomAPI_Interpolate builder(_handle, periodic, tol); +// cout << "done" << endl; +// auto curve = builder.Curve(); +// cout << "curve done" << endl; +// return BRepBuilderAPI_MakeEdge(curve); +//// return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); +// }, +// py::arg("points"), +// py::arg("periodic")=false, +// py::arg("tol")=1e-8, +// R"delimiter( +//Generate a piecewise continuous spline-curve interpolating a list of points. +// +//Parameters +//---------- +// +//points : List[gp_Pnt] or Tuple[gp_Pnt] or np.ndarray[double] +// List (or tuple) of gp_Pnt. If a numpy array is provided instead, the data must contain the coordinates +// +//periodic : bool +// Whether the result should be periodic +// +//tol : float +// Tolerance for the distance between points. +// +//)delimiter"); + + + m.def("SplineSurfaceApproximation", [](py::array_t pnt_array, + Approx_ParametrizationType approx_type, int deg_min, int deg_max, GeomAbs_Shape continuity, double tol, + bool periodic, double degen_tol) { + if (pnt_array.ndim() != 3) + throw Exception("`points` array must have dimension 3."); + if (pnt_array.shape(2) != 3) + throw Exception("The third dimension must have size 3."); + + auto array = py::extract>(pnt_array)(); + TColgp_Array2OfPnt points(1, pnt_array.shape(0), 1, pnt_array.shape(1)); + auto pnts_unchecked = pnt_array.unchecked<3>(); + for (int i = 0; i < pnt_array.shape(0); ++i) + for (int j = 0; j < pnt_array.shape(1); ++j) + points.SetValue(i+1, j+1, gp_Pnt(pnts_unchecked(i, j, 0), pnts_unchecked(i, j, 1), pnts_unchecked(i, j, 2))); + + GeomAPI_PointsToBSplineSurface builder; + builder.Init(points, approx_type, deg_min, deg_max, continuity, tol, periodic); + return BRepBuilderAPI_MakeFace(builder.Surface(), tol).Face(); + }, + py::arg("points"), + py::arg("approx_type") = Approx_ParametrizationType::Approx_ChordLength, + py::arg("deg_min") = 3, + py::arg("deg_max") = 8, + py::arg("continuity") = GeomAbs_Shape::GeomAbs_C2, + py::arg("tol") = 1e-3, + py::arg("periodic") = false, + py::arg("degen_tol") = 1e-8, + R"delimiter( +Generate a piecewise continuous spline-surface approximating an array of points. + +Parameters +---------- + +points : np.ndarray + Array of points coordinates. The first dimension corresponds to the first surface coordinate point + index, the second dimension to the second surface coordinate point index. The third dimension refers to physical + coordinates. Such an array can be generated with code like:: + + px, py = np.meshgrid(*[np.linspace(0, 1, N)]*2) + points = np.array([[(px[i,j], py[i,j], px[i,j]*py[i,j]**2) for j in range(N)] for i in range(N)]) + +approx_type : ApproxParamType + Assumption on location of parameters wrt points. + +deg_min : int + Minimum polynomial degree of splines + +deg_max : int + Maxmium polynomial degree of splines + +continuity : ShapeContinuity + Continuity requirement on the approximating surface + +tol : float + Tolerance for the distance from individual points to the approximating surface. + +periodic : bool + Whether the result should be periodic + +degen_tol : double + Tolerance for resolution of degenerate edges + +)delimiter"); + + m.def("SplineSurfaceInterpolation", []( + py::array_t pnt_array, Approx_ParametrizationType approx_type, bool periodic, double degen_tol) { + + if (pnt_array.ndim() != 3) + throw Exception("`points` array must have dimension 3."); + if (pnt_array.shape(2) != 3) + throw Exception("The third dimension must have size 3."); + + auto array = py::extract>(pnt_array)(); + TColgp_Array2OfPnt points(1, pnt_array.shape(0), 1, pnt_array.shape(1)); + auto pnts_unchecked = pnt_array.unchecked<3>(); + for (int i = 0; i < pnt_array.shape(0); ++i) + for (int j = 0; j < pnt_array.shape(1); ++j) + points.SetValue(i+1, j+1, gp_Pnt(pnts_unchecked(i, j, 0), pnts_unchecked(i, j, 1), pnts_unchecked(i, j, 2))); + + GeomAPI_PointsToBSplineSurface builder; + builder.Interpolate(points, approx_type, periodic); + return BRepBuilderAPI_MakeFace(builder.Surface(), degen_tol).Face(); + }, + py::arg("points"), + py::arg("approx_type") = Approx_ParametrizationType::Approx_ChordLength, + py::arg("periodic") = false, + py::arg("degen_tol") = 1e-8, + R"delimiter( +Generate a piecewise continuous spline-surface interpolating an array of points. + +Parameters +---------- + +points : np.ndarray + Array of points coordinates. The first dimension corresponds to the first surface coordinate point + index, the second dimension to the second surface coordinate point index. The third dimension refers to physical + coordinates. Such an array can be generated with code like:: + + px, py = np.meshgrid(*[np.linspace(0, 1, N)]*2) + points = np.array([[(px[i,j], py[i,j], px[i,j]*py[i,j]**2) for j in range(N)] for i in range(N)]) + +approx_type : ApproxParamType + Assumption on location of parameters wrt points. + +periodic : bool + Whether the result should be periodic + +degen_tol : double + Tolerance for resolution of degenerate edges + +)delimiter"); m.def("MakeFillet", [](TopoDS_Shape shape, std::vector edges, double r) { From a5aed39f9dd4888bc1c4142d84e22a334650e35f Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Thu, 2 Dec 2021 16:52:38 +0100 Subject: [PATCH 02/20] SplineInterpolation now works; refined some docstrings --- libsrc/occ/python_occ_shapes.cpp | 101 ++++++++++++++----------------- 1 file changed, 45 insertions(+), 56 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index f610fdbf..d9803f20 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -2069,60 +2069,49 @@ tol : float )delimiter"); -// CRASHES BADLY during call of Curve()! -// m.def("SplineInterpolation", [](py::object pnts, bool periodic, double tol) { -// cout << "enter" << endl; -// TColgp_Array1OfPnt points(0, 0); -// cout << "points array exists" << endl; -// if (py::extract>(pnts).check()) { -// std::vector pnt_list{py::extract>(pnts)()}; -// points.Resize(1, pnt_list.size(), true); -// for (int i = 0; i < pnt_list.size(); i++) -// points.SetValue(i+1, pnt_list[i]); -// } else if (py::extract>(pnts).check()) { -// py::array_t pnt_array{py::extract>(pnts)()}; -// if (pnt_array.ndim() != 2) -// throw Exception("`points` array must have dimension 2."); -// if (pnt_array.shape(1) != 3) -// throw Exception("The second dimension must have size 3."); -// cout << "resize" << endl; -// points.Resize(1, pnt_array.shape(0), true); -// auto pnts_unchecked = pnt_array.unchecked<2>(); -// for (int i = 0; i < pnt_array.shape(0); ++i) -// points.SetValue(i+1, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); -// cout << "values set" << endl; -// } else -// throw Exception("Not able to process the data type of points"); -// -// TColgp_HArray1OfPnt hpoints{points}; -// const auto _handle = opencascade::handle{&hpoints}; -// cout << "build" << endl; -// GeomAPI_Interpolate builder(_handle, periodic, tol); -// cout << "done" << endl; -// auto curve = builder.Curve(); -// cout << "curve done" << endl; -// return BRepBuilderAPI_MakeEdge(curve); -//// return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); -// }, -// py::arg("points"), -// py::arg("periodic")=false, -// py::arg("tol")=1e-8, -// R"delimiter( -//Generate a piecewise continuous spline-curve interpolating a list of points. -// -//Parameters -//---------- -// -//points : List[gp_Pnt] or Tuple[gp_Pnt] or np.ndarray[double] -// List (or tuple) of gp_Pnt. If a numpy array is provided instead, the data must contain the coordinates -// -//periodic : bool -// Whether the result should be periodic -// -//tol : float -// Tolerance for the distance between points. -// -//)delimiter"); + m.def("SplineInterpolation", [](py::object pnts, bool periodic, double tol) { + Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, 1); + if (py::extract>(pnts).check()) { + std::vector pnt_list{py::extract>(pnts)()}; + points->Resize(1, pnt_list.size(), true); + for (int i = 0; i < pnt_list.size(); i++) + points->SetValue(i+1, pnt_list[i]); + } else if (py::extract>(pnts).check()) { + py::array_t pnt_array{py::extract>(pnts)()}; + if (pnt_array.ndim() != 2) + throw Exception("`points` array must have dimension 2."); + if (pnt_array.shape(1) != 3) + throw Exception("The second dimension must have size 3."); + points->Resize(1, pnt_array.shape(0), false); + auto pnts_unchecked = pnt_array.unchecked<2>(); + for (int i = 0; i < pnt_array.shape(0); ++i) + points->SetValue(i+1, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); + } else + throw Exception("Not able to process the data type of points"); + + GeomAPI_Interpolate builder(points, periodic, tol); + builder.Perform(); + return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); + }, + py::arg("points"), + py::arg("periodic")=false, + py::arg("tol")=1e-8, + R"delimiter( +Generate a piecewise continuous spline-curve interpolating a list of points. + +Parameters +---------- + +points : List[gp_Pnt] or Tuple[gp_Pnt] or np.ndarray[double] + List (or tuple) of gp_Pnt. If a numpy array is provided instead, the data must contain the coordinates + +periodic : bool + Whether the result should be periodic + +tol : float + Tolerance for the distance between points. + +)delimiter"); m.def("SplineSurfaceApproximation", [](py::array_t pnt_array, @@ -2182,7 +2171,7 @@ tol : float Tolerance for the distance from individual points to the approximating surface. periodic : bool - Whether the result should be periodic + Whether the result should be periodic in the first surface parameter degen_tol : double Tolerance for resolution of degenerate edges @@ -2230,7 +2219,7 @@ approx_type : ApproxParamType Assumption on location of parameters wrt points. periodic : bool - Whether the result should be periodic + Whether the result should be periodic in the first surface parameter degen_tol : double Tolerance for resolution of degenerate edges From 193a7001e44c925170b38e2436a389cd73570062 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Thu, 2 Dec 2021 17:39:11 +0100 Subject: [PATCH 03/20] added 2d spline interpolation but not added to workplane yet --- libsrc/occ/python_occ_shapes.cpp | 73 ++++++++++++++++++++++++-------- 1 file changed, 55 insertions(+), 18 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index d9803f20..c82f2e90 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -33,16 +33,19 @@ #include #include #include +#include #include #include #include #include #include #include +#include #include #include #include #include +#include #include #include #include @@ -2012,6 +2015,7 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return BRepBuilderAPI_MakeEdge(curve).Edge(); }, py::arg("points"), "create Bezier curve"); + //TODO: 2d version m.def("SplineApproximation", [](py::object pnts, Approx_ParametrizationType approx_type, int deg_min, int deg_max, GeomAbs_Shape continuity, double tol) { TColgp_Array1OfPnt points(0, 0); @@ -2070,28 +2074,61 @@ tol : float )delimiter"); m.def("SplineInterpolation", [](py::object pnts, bool periodic, double tol) { - Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, 1); - if (py::extract>(pnts).check()) { + + auto _2d = [](const Handle(TColgp_HArray1OfPnt2d) &points, bool periodic, double tol) { + Geom2dAPI_Interpolate builder(points, periodic, tol); + builder.Perform(); + return BRepBuilderAPI_MakeEdge2d(builder.Curve()).Edge(); + }; + + auto _3d = [](const Handle(TColgp_HArray1OfPnt) &points, bool periodic, double tol) { + GeomAPI_Interpolate builder(points, periodic, tol); + builder.Perform(); + return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); + }; + + if (py::extract>(pnts).check()) + { std::vector pnt_list{py::extract>(pnts)()}; - points->Resize(1, pnt_list.size(), true); + Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt_list.size()); for (int i = 0; i < pnt_list.size(); i++) points->SetValue(i+1, pnt_list[i]); - } else if (py::extract>(pnts).check()) { + return _3d(points, periodic, tol); + } + else if(py::extract>(pnts).check()) + { + std::vector pnt_list{py::extract>(pnts)()}; + Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d(1, pnt_list.size()); + for (int i = 0; i < pnt_list.size(); i++) + points->SetValue(i+1, pnt_list[i]); + return _2d(points, periodic, tol); + } + else if (py::extract>(pnts).check()) + { py::array_t pnt_array{py::extract>(pnts)()}; if (pnt_array.ndim() != 2) throw Exception("`points` array must have dimension 2."); - if (pnt_array.shape(1) != 3) - throw Exception("The second dimension must have size 3."); - points->Resize(1, pnt_array.shape(0), false); - auto pnts_unchecked = pnt_array.unchecked<2>(); - for (int i = 0; i < pnt_array.shape(0); ++i) - points->SetValue(i+1, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); - } else - throw Exception("Not able to process the data type of points"); - - GeomAPI_Interpolate builder(points, periodic, tol); - builder.Perform(); - return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); + if (pnt_array.shape(1) == 3) + { + Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt_array.shape(0)); + auto pnts_unchecked = pnt_array.unchecked<2>(); + for (int i = 0; i < pnt_array.shape(0); ++i) + points->SetValue(i+1, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); + return _3d(points, periodic, tol); + } + else if (pnt_array.shape(1) == 2) + { + Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d(1, pnt_array.shape(0)); + auto pnts_unchecked = pnt_array.unchecked<2>(); + for (int i = 0; i < pnt_array.shape(0); ++i) + points->SetValue(i+1, gp_Pnt2d(pnts_unchecked(i, 0), pnts_unchecked(i, 1))); + return _2d(points, periodic, tol); + } + else + throw Exception("The second dimension must have size 2 or 3, but has " + to_string(pnt_array.shape(1))); + } + else + throw Exception("Not able to process the data type of points"); }, py::arg("points"), py::arg("periodic")=false, @@ -2102,8 +2139,8 @@ Generate a piecewise continuous spline-curve interpolating a list of points. Parameters ---------- -points : List[gp_Pnt] or Tuple[gp_Pnt] or np.ndarray[double] - List (or tuple) of gp_Pnt. If a numpy array is provided instead, the data must contain the coordinates +points : List|Tuple[gp_Pnt|gp_Pnt2d] or np.ndarray[double] + List (or tuple) of gp_Pnt (or gp_Pnt2d). If a numpy array is provided instead, the data must contain the coordinates periodic : bool Whether the result should be periodic From 31fa22626cf13f08363246bd9f35d3a5b1c197a7 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Fri, 3 Dec 2021 11:55:02 +0100 Subject: [PATCH 04/20] Split 2d and 3d spline implementations, use tangent data; added Spline member to WorkPlane --- libsrc/occ/python_occ_shapes.cpp | 284 ++++++++++++++++++++++--------- 1 file changed, 205 insertions(+), 79 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index c82f2e90..fe5a0ecd 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include @@ -442,6 +443,75 @@ public: return shared_from_this(); } + auto Spline(const std::vector &points, bool periodic, double tol, const std::map &tangents, optional name = nullopt) + { + gp_Pnt2d P1 = localpos.Location(); + gp_Pnt P13d = surf->Value(P1.X(), P1.Y()); + + gp_Pnt2d PLast = points.back(); + gp_Pnt PLast3d = surf->Value(PLast.X(), PLast.Y()); + + Handle(TColgp_HArray1OfPnt2d) allpoints = new TColgp_HArray1OfPnt2d(1, points.size() + 1); + allpoints->SetValue(1, P1); + for (int i = 0; i < points.size(); i++) + allpoints->SetValue(i+2, points[i]); + Geom2dAPI_Interpolate builder(allpoints, periodic, tol); + + if (tangents.size() > 0) + { + const gp_Vec2d dummy_vec = tangents.begin()->second; + TColgp_Array1OfVec2d tangent_vecs(1, allpoints->Length()); + Handle(TColStd_HArray1OfBoolean) tangent_flags = new TColStd_HArray1OfBoolean(1, allpoints->Length()); + for (int i : Range(allpoints->Length())) + { + if (tangents.count(i) > 0) + { + tangent_vecs.SetValue(i+1, tangents.at(i)); + tangent_flags->SetValue(i+1, true); + } + else + { + tangent_vecs.SetValue(i+1, dummy_vec); + tangent_flags->SetValue(i+1, false); + } + } + builder.Load(tangent_vecs, tangent_flags); + } + + + builder.Perform(); + auto curve2d = builder.Curve(); + + if (startvertex.IsNull()) + startvertex = lastvertex = BRepBuilderAPI_MakeVertex(P13d).Vertex(); + auto endv = periodic ? startvertex : BRepBuilderAPI_MakeVertex(PLast3d).Vertex(); + + //create 3d edge from 2d curve using surf + auto edge = BRepBuilderAPI_MakeEdge(curve2d, surf, lastvertex, endv).Edge(); + lastvertex = endv; + BRepLib::BuildCurves3d(edge); + wire_builder.Add(edge); + + // update localpos + localpos.SetLocation(PLast); + + //compute angle of rotation + //compute tangent t2 in PLast + const auto dir = localpos.Direction(); + gp_Vec2d t = gp_Vec2d(dir.X(), dir.Y()); + gp_Vec2d t2 = curve2d->DN(curve2d->LastParameter(), 1); + + double angle = t.Angle(t2); //angle \in [-pi,pi] + + //update localpos.Direction() + Rotate(angle*180/M_PI); + + if (periodic) + Finish(); + + return shared_from_this(); + } + auto ArcTo (double h, double v, const gp_Vec2d t) { gp_Pnt2d P1 = localpos.Location(); @@ -1855,6 +1925,100 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return curve; */ }, py::arg("c"), py::arg("r"), "create 2d circle curve"); + + m.def("SplineApproximation", [](const std::vector &points, Approx_ParametrizationType approx_type, int deg_min, + int deg_max, GeomAbs_Shape continuity, double tol) -> Handle(Geom2d_Curve) { + TColgp_Array1OfPnt2d hpoints(0, 0); + hpoints.Resize(0, points.size() - 1, true); + for (int i = 0; i < points.size(); i++) + hpoints.SetValue(i, points[i]); + + Geom2dAPI_PointsToBSpline builder(hpoints, approx_type, deg_min, deg_max, continuity, tol); + return Handle(Geom2d_BSplineCurve)(builder.Curve()); + }, + py::arg("points"), + py::arg("approx_type") = Approx_ParametrizationType::Approx_ChordLength, + py::arg("deg_min") = 3, + py::arg("deg_max") = 8, + py::arg("continuity") = GeomAbs_Shape::GeomAbs_C2, + py::arg("tol")=1e-8, + R"delimiter( +Generate a piecewise continuous spline-curve approximating a list of points in 2d. + +Parameters +---------- + +points : List|Tuple[gp_Pnt2d] + List (or tuple) of gp_Pnt. + +approx_type : ApproxParamType + Assumption on location of parameters wrt points. + +deg_min : int + Minimum polynomial degree of splines + +deg_max : int + Maxmium polynomial degree of splines + +continuity : ShapeContinuity + Continuity requirement on the approximating surface + +tol : float + Tolerance for the distance from individual points to the approximating curve. + +)delimiter"); + + m.def("SplineInterpolation", [](const std::vector &points, bool periodic, double tol, const std::map &tangents) -> Handle(Geom2d_Curve) { + Handle(TColgp_HArray1OfPnt2d) hpoints = new TColgp_HArray1OfPnt2d(1, points.size()); + for (int i = 0; i < points.size(); i++) + hpoints->SetValue(i+1, points[i]); + Geom2dAPI_Interpolate builder(hpoints, periodic, tol); + + if (tangents.size() > 0) + { + const gp_Vec2d dummy_vec = tangents.begin()->second; + TColgp_Array1OfVec2d tangent_vecs(1, points.size()); + Handle(TColStd_HArray1OfBoolean) tangent_flags = new TColStd_HArray1OfBoolean(1, points.size()); + for (int i : Range(points.size())) + { + if (tangents.count(i) > 0) + { + tangent_vecs.SetValue(i+1, tangents.at(i)); + tangent_flags->SetValue(i+1, true); + } else{ + tangent_vecs.SetValue(i+1, dummy_vec); + tangent_flags->SetValue(i+1, false); + } + } + builder.Load(tangent_vecs, tangent_flags); + } + + builder.Perform(); + return Handle(Geom2d_BSplineCurve)(builder.Curve()); + }, + py::arg("points"), + py::arg("periodic")=false, + py::arg("tol")=1e-8, + py::arg("tangents")=std::map{}, + R"delimiter( +Generate a piecewise continuous spline-curve interpolating a list of points in 2d. + +Parameters +---------- + +points : List|Tuple[gp_Pnt2d] + List (or tuple) of gp_Pnt2d. + +periodic : bool + Whether the result should be periodic + +tol : float + Tolerance for the distance between points. + +tangents : Dict[int, gp_Vec2d] + Tangent vectors for the points indicated by the key value (0-based). + +)delimiter"); m.def("Glue", [] (const std::vector shapes) -> TopoDS_Shape { @@ -2015,30 +2179,15 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) return BRepBuilderAPI_MakeEdge(curve).Edge(); }, py::arg("points"), "create Bezier curve"); - //TODO: 2d version - m.def("SplineApproximation", [](py::object pnts, Approx_ParametrizationType approx_type, int deg_min, + + m.def("SplineApproximation", [](const std::vector &points, Approx_ParametrizationType approx_type, int deg_min, int deg_max, GeomAbs_Shape continuity, double tol) { - TColgp_Array1OfPnt points(0, 0); - if (py::extract>(pnts).check()) { - std::vector pnt_list{py::extract>(pnts)()}; - points.Resize(0, pnt_list.size()-1, true); - for (int i = 0; i < pnt_list.size(); i++) - points.SetValue(i, pnt_list[i]); - } else if (py::extract>(pnts).check()) { - py::array_t pnt_array{py::extract>(pnts)()}; - if (pnt_array.ndim() != 2) - throw Exception("`points` array must have dimension 2."); - if (pnt_array.shape(1) != 3) - throw Exception("The second dimension must have size 3."); + TColgp_Array1OfPnt hpoints(0, 0); + hpoints.Resize(0, points.size() - 1, true); + for (int i = 0; i < points.size(); i++) + hpoints.SetValue(i, points[i]); - points.Resize(0, pnt_array.shape(0)-1, true); - auto pnts_unchecked = pnt_array.unchecked<2>(); - for (int i = 0; i < pnt_array.shape(0); ++i) - points.SetValue(i, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); - } else - throw Exception("Not able to process the data type of points"); - - GeomAPI_PointsToBSpline builder(points, approx_type, deg_min, deg_max, continuity, tol); + GeomAPI_PointsToBSpline builder(hpoints, approx_type, deg_min, deg_max, continuity, tol); return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); }, py::arg("points"), @@ -2048,13 +2197,13 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) py::arg("continuity") = GeomAbs_Shape::GeomAbs_C2, py::arg("tol")=1e-8, R"delimiter( -Generate a piecewise continuous spline-curve approximating a list of points. +Generate a piecewise continuous spline-curve approximating a list of points in 3d. Parameters ---------- -points : List[gp_Pnt] or Tuple[gp_Pnt] or np.ndarray[double] - List (or tuple) of gp_Pnt. If a numpy array is provided instead, the data must contain the coordinates +points : List[gp_Pnt] or Tuple[gp_Pnt] + List (or tuple) of gp_Pnt. approx_type : ApproxParamType Assumption on location of parameters wrt points. @@ -2073,74 +2222,47 @@ tol : float )delimiter"); - m.def("SplineInterpolation", [](py::object pnts, bool periodic, double tol) { + m.def("SplineInterpolation", [](const std::vector &points, bool periodic, double tol, const std::map &tangents) { + Handle(TColgp_HArray1OfPnt) hpoints = new TColgp_HArray1OfPnt(1, points.size()); + for (int i = 0; i < points.size(); i++) + hpoints->SetValue(i+1, points[i]); - auto _2d = [](const Handle(TColgp_HArray1OfPnt2d) &points, bool periodic, double tol) { - Geom2dAPI_Interpolate builder(points, periodic, tol); - builder.Perform(); - return BRepBuilderAPI_MakeEdge2d(builder.Curve()).Edge(); - }; + GeomAPI_Interpolate builder(hpoints, periodic, tol); - auto _3d = [](const Handle(TColgp_HArray1OfPnt) &points, bool periodic, double tol) { - GeomAPI_Interpolate builder(points, periodic, tol); - builder.Perform(); - return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); - }; - - if (py::extract>(pnts).check()) + if (tangents.size() > 0) + { + const gp_Vec dummy_vec = tangents.begin()->second; + TColgp_Array1OfVec tangent_vecs(1, points.size()); + Handle(TColStd_HArray1OfBoolean) tangent_flags = new TColStd_HArray1OfBoolean(1, points.size()); + for (int i : Range(points.size())) { - std::vector pnt_list{py::extract>(pnts)()}; - Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt_list.size()); - for (int i = 0; i < pnt_list.size(); i++) - points->SetValue(i+1, pnt_list[i]); - return _3d(points, periodic, tol); + if (tangents.count(i) > 0) + { + tangent_vecs.SetValue(i+1, tangents.at(i)); + tangent_flags->SetValue(i+1, true); + } else{ + tangent_vecs.SetValue(i+1, dummy_vec); + tangent_flags->SetValue(i+1, false); + } } - else if(py::extract>(pnts).check()) - { - std::vector pnt_list{py::extract>(pnts)()}; - Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d(1, pnt_list.size()); - for (int i = 0; i < pnt_list.size(); i++) - points->SetValue(i+1, pnt_list[i]); - return _2d(points, periodic, tol); - } - else if (py::extract>(pnts).check()) - { - py::array_t pnt_array{py::extract>(pnts)()}; - if (pnt_array.ndim() != 2) - throw Exception("`points` array must have dimension 2."); - if (pnt_array.shape(1) == 3) - { - Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, pnt_array.shape(0)); - auto pnts_unchecked = pnt_array.unchecked<2>(); - for (int i = 0; i < pnt_array.shape(0); ++i) - points->SetValue(i+1, gp_Pnt(pnts_unchecked(i, 0), pnts_unchecked(i, 1), pnts_unchecked(i, 2))); - return _3d(points, periodic, tol); - } - else if (pnt_array.shape(1) == 2) - { - Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d(1, pnt_array.shape(0)); - auto pnts_unchecked = pnt_array.unchecked<2>(); - for (int i = 0; i < pnt_array.shape(0); ++i) - points->SetValue(i+1, gp_Pnt2d(pnts_unchecked(i, 0), pnts_unchecked(i, 1))); - return _2d(points, periodic, tol); - } - else - throw Exception("The second dimension must have size 2 or 3, but has " + to_string(pnt_array.shape(1))); + builder.Load(tangent_vecs, tangent_flags); } - else - throw Exception("Not able to process the data type of points"); + + builder.Perform(); + return BRepBuilderAPI_MakeEdge(builder.Curve()).Edge(); }, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, + py::arg("tangents")=std::map{}, R"delimiter( -Generate a piecewise continuous spline-curve interpolating a list of points. +Generate a piecewise continuous spline-curve interpolating a list of points in 3d. Parameters ---------- -points : List|Tuple[gp_Pnt|gp_Pnt2d] or np.ndarray[double] - List (or tuple) of gp_Pnt (or gp_Pnt2d). If a numpy array is provided instead, the data must contain the coordinates +points : List|Tuple[gp_Pnt] + List (or tuple) of gp_Pnt periodic : bool Whether the result should be periodic @@ -2148,6 +2270,9 @@ periodic : bool tol : float Tolerance for the distance between points. +tangents : Dict[int, gp_Vec] + Tangent vectors for the points indicated by the key value (0-based). + )delimiter"); @@ -2324,6 +2449,7 @@ degen_tol : double py::arg("l"), py::arg("name")=nullopt) .def("Line", [](WorkPlane&wp,double h,double v, optional name) { return wp.Line(h,v,name); }, py::arg("dx"), py::arg("dy"), py::arg("name")=nullopt) + .def("Spline", &WorkPlane::Spline, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, py::arg("tangents")=std::map{}, py::arg("name")=nullopt, "draw spline starting from current position, tangents can be given for each point (0 refers to current position)") .def("Rectangle", &WorkPlane::Rectangle, py::arg("l"), py::arg("w"), "draw rectangle, with current position as corner, use current direction") .def("RectangleC", &WorkPlane::RectangleCentered, py::arg("l"), py::arg("w"), "draw rectangle, with current position as center, use current direction") .def("Circle", [](WorkPlane&wp, double x, double y, double r) { From 9f83730fb5ae8f7e2e2ea947da3e32a373ba18b1 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Mon, 6 Dec 2021 15:26:57 +0100 Subject: [PATCH 05/20] add a check on first point given to WP::Spline; more precise docs --- libsrc/occ/python_occ_shapes.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index fe5a0ecd..28983a71 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -451,10 +451,13 @@ public: gp_Pnt2d PLast = points.back(); gp_Pnt PLast3d = surf->Value(PLast.X(), PLast.Y()); + if (points.front().Distance(P1) <= tol) + throw Exception("First item of given list of points is too close to current position (distance <= tol).") + Handle(TColgp_HArray1OfPnt2d) allpoints = new TColgp_HArray1OfPnt2d(1, points.size() + 1); allpoints->SetValue(1, P1); for (int i = 0; i < points.size(); i++) - allpoints->SetValue(i+2, points[i]); + allpoints->SetValue(i + 2, points[i]); Geom2dAPI_Interpolate builder(allpoints, periodic, tol); if (tangents.size() > 0) @@ -2449,7 +2452,9 @@ degen_tol : double py::arg("l"), py::arg("name")=nullopt) .def("Line", [](WorkPlane&wp,double h,double v, optional name) { return wp.Line(h,v,name); }, py::arg("dx"), py::arg("dy"), py::arg("name")=nullopt) - .def("Spline", &WorkPlane::Spline, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, py::arg("tangents")=std::map{}, py::arg("name")=nullopt, "draw spline starting from current position, tangents can be given for each point (0 refers to current position)") + .def("Spline", &WorkPlane::Spline, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, + py::arg("tangents")=std::map{}, py::arg("name")=nullopt, + "draw spline starting from current position (implicitly added to given list of points), tangents can be specified for each point (0 refers to current position)") .def("Rectangle", &WorkPlane::Rectangle, py::arg("l"), py::arg("w"), "draw rectangle, with current position as corner, use current direction") .def("RectangleC", &WorkPlane::RectangleCentered, py::arg("l"), py::arg("w"), "draw rectangle, with current position as center, use current direction") .def("Circle", [](WorkPlane&wp, double x, double y, double r) { From fdf26641dde3d93fbf1bda7315e90aa8f67d3835 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Mon, 6 Dec 2021 16:28:02 +0100 Subject: [PATCH 06/20] fixed exception --- libsrc/occ/python_occ_shapes.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 28983a71..f4fd387f 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -452,7 +452,7 @@ public: gp_Pnt PLast3d = surf->Value(PLast.X(), PLast.Y()); if (points.front().Distance(P1) <= tol) - throw Exception("First item of given list of points is too close to current position (distance <= tol).") + throw Exception("First item of given list of points is too close to current position (distance <= tol)."); Handle(TColgp_HArray1OfPnt2d) allpoints = new TColgp_HArray1OfPnt2d(1, points.size() + 1); allpoints->SetValue(1, P1); From 6656181e2b5ba2892e58ca51d727c1d23d243dfb Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Wed, 22 Dec 2021 16:21:24 +0100 Subject: [PATCH 07/20] work om WP Spline: detect closing similar to ArcTo and LineTo, remove "name" arg --- libsrc/occ/python_occ_shapes.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index f4fd387f..6e11fffb 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -443,7 +443,7 @@ public: return shared_from_this(); } - auto Spline(const std::vector &points, bool periodic, double tol, const std::map &tangents, optional name = nullopt) + auto Spline(const std::vector &points, bool periodic, double tol, const std::map &tangents) { gp_Pnt2d P1 = localpos.Location(); gp_Pnt P13d = surf->Value(P1.X(), P1.Y()); @@ -485,9 +485,10 @@ public: builder.Perform(); auto curve2d = builder.Curve(); + const bool closing = periodic || PLast.Distance(startpnt) < 1e-10; if (startvertex.IsNull()) startvertex = lastvertex = BRepBuilderAPI_MakeVertex(P13d).Vertex(); - auto endv = periodic ? startvertex : BRepBuilderAPI_MakeVertex(PLast3d).Vertex(); + auto endv = closing ? startvertex : BRepBuilderAPI_MakeVertex(PLast3d).Vertex(); //create 3d edge from 2d curve using surf auto edge = BRepBuilderAPI_MakeEdge(curve2d, surf, lastvertex, endv).Edge(); @@ -509,7 +510,7 @@ public: //update localpos.Direction() Rotate(angle*180/M_PI); - if (periodic) + if (closing) Finish(); return shared_from_this(); @@ -2453,7 +2454,7 @@ degen_tol : double .def("Line", [](WorkPlane&wp,double h,double v, optional name) { return wp.Line(h,v,name); }, py::arg("dx"), py::arg("dy"), py::arg("name")=nullopt) .def("Spline", &WorkPlane::Spline, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, - py::arg("tangents")=std::map{}, py::arg("name")=nullopt, + py::arg("tangents")=std::map{}, "draw spline starting from current position (implicitly added to given list of points), tangents can be specified for each point (0 refers to current position)") .def("Rectangle", &WorkPlane::Rectangle, py::arg("l"), py::arg("w"), "draw rectangle, with current position as corner, use current direction") .def("RectangleC", &WorkPlane::RectangleCentered, py::arg("l"), py::arg("w"), "draw rectangle, with current position as center, use current direction") From 88bb6ec6af947b1a0e4c77f03e2e517422f4a192 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Mon, 27 Dec 2021 16:10:37 +0100 Subject: [PATCH 08/20] add some functions to access WorkPlane data and the possibility to create splines from any starting point --- libsrc/occ/python_occ_shapes.cpp | 46 +++++++++++++++++++++++++------- 1 file changed, 37 insertions(+), 9 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 6e11fffb..eb26d1ed 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -362,6 +362,20 @@ public: return shared_from_this(); } + auto StartPnt() const { + return startpnt; + } + + auto CurrentLocation() const + { + return localpos.Location(); + } + + auto CurrentDirection() const + { + return gp_Vec2d(localpos.Direction()); + } + auto MoveTo (double h, double v) { startpnt = gp_Pnt2d(h,v); @@ -443,21 +457,32 @@ public: return shared_from_this(); } - auto Spline(const std::vector &points, bool periodic, double tol, const std::map &tangents) + auto Spline(const std::vector &points, bool periodic, double tol, const std::map &tangents, + bool start_from_localpos) { - gp_Pnt2d P1 = localpos.Location(); + gp_Pnt2d P1 = start_from_localpos ? localpos.Location() : points.front(); gp_Pnt P13d = surf->Value(P1.X(), P1.Y()); gp_Pnt2d PLast = points.back(); gp_Pnt PLast3d = surf->Value(PLast.X(), PLast.Y()); - if (points.front().Distance(P1) <= tol) - throw Exception("First item of given list of points is too close to current position (distance <= tol)."); + Handle(TColgp_HArray1OfPnt2d) allpoints; + if (start_from_localpos) + { + if (points.front().Distance(P1) <= tol) + throw Exception("First item of given list of points is too close to current position (distance <= tol)."); + allpoints = new TColgp_HArray1OfPnt2d(1, points.size() + 1); + allpoints->SetValue(1, P1); + for (int i = 0; i < points.size(); i++) + allpoints->SetValue(i + 2, points[i]); + } + else + { + allpoints = new TColgp_HArray1OfPnt2d(1, points.size()); + for (int i = 0; i < points.size(); i++) + allpoints->SetValue(i + 1, points[i]); + } - Handle(TColgp_HArray1OfPnt2d) allpoints = new TColgp_HArray1OfPnt2d(1, points.size() + 1); - allpoints->SetValue(1, P1); - for (int i = 0; i < points.size(); i++) - allpoints->SetValue(i + 2, points[i]); Geom2dAPI_Interpolate builder(allpoints, periodic, tol); if (tangents.size() > 0) @@ -2440,6 +2465,9 @@ degen_tol : double py::class_> (m, "WorkPlane") .def(py::init(), py::arg("axes")=gp_Ax3(), py::arg("pos")=gp_Ax2d()) + .def_property_readonly("CurrentLocation", &WorkPlane::CurrentLocation) + .def_property_readonly("CurrentDirection", &WorkPlane::CurrentDirection) + .def_property_readonly("StartPnt", &WorkPlane::StartPnt) .def("MoveTo", &WorkPlane::MoveTo, py::arg("h"), py::arg("v"), "moveto (h,v), and start new wire") .def("Move", &WorkPlane::Move, py::arg("l"), "move 'l' from current position and direction, start new wire") .def("Direction", &WorkPlane::Direction, py::arg("dirh"), py::arg("dirv"), "reset direction to (dirh, dirv)") @@ -2454,7 +2482,7 @@ degen_tol : double .def("Line", [](WorkPlane&wp,double h,double v, optional name) { return wp.Line(h,v,name); }, py::arg("dx"), py::arg("dy"), py::arg("name")=nullopt) .def("Spline", &WorkPlane::Spline, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, - py::arg("tangents")=std::map{}, + py::arg("tangents")=std::map{}, py::arg("start_from_localpos")=true, "draw spline starting from current position (implicitly added to given list of points), tangents can be specified for each point (0 refers to current position)") .def("Rectangle", &WorkPlane::Rectangle, py::arg("l"), py::arg("w"), "draw rectangle, with current position as corner, use current direction") .def("RectangleC", &WorkPlane::RectangleCentered, py::arg("l"), py::arg("w"), "draw rectangle, with current position as center, use current direction") From 093692f825721a47fcded672dcde8bf9ec1465d5 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Mon, 27 Dec 2021 17:08:00 +0100 Subject: [PATCH 09/20] docstring fix --- libsrc/occ/python_occ_shapes.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index eb26d1ed..126aa23e 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -2483,7 +2483,7 @@ degen_tol : double py::arg("dx"), py::arg("dy"), py::arg("name")=nullopt) .def("Spline", &WorkPlane::Spline, py::arg("points"), py::arg("periodic")=false, py::arg("tol")=1e-8, py::arg("tangents")=std::map{}, py::arg("start_from_localpos")=true, - "draw spline starting from current position (implicitly added to given list of points), tangents can be specified for each point (0 refers to current position)") + "draw spline (default: starting from current position, which is implicitly added to given list of points), tangents can be specified for each point (0 refers to starting point)") .def("Rectangle", &WorkPlane::Rectangle, py::arg("l"), py::arg("w"), "draw rectangle, with current position as corner, use current direction") .def("RectangleC", &WorkPlane::RectangleCentered, py::arg("l"), py::arg("w"), "draw rectangle, with current position as center, use current direction") .def("Circle", [](WorkPlane&wp, double x, double y, double r) { From 51013105edd572152346a4f0aa9d649d4ed96282 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Tue, 8 Feb 2022 14:01:25 +0100 Subject: [PATCH 10/20] doc typo --- libsrc/occ/python_occ_shapes.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 679c3bb7..282044c6 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -1870,7 +1870,7 @@ deg_min : int Minimum polynomial degree of splines deg_max : int - Maxmium polynomial degree of splines + Maximum polynomial degree of splines continuity : ShapeContinuity Continuity requirement on the approximating surface From 5f6b23426231199c1bb016e8d5a6101e84185e15 Mon Sep 17 00:00:00 2001 From: Matthias Rambausek Date: Fri, 11 Feb 2022 11:46:55 +0100 Subject: [PATCH 11/20] shorter property names --- libsrc/occ/python_occ_shapes.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 282044c6..57479a4b 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -2348,9 +2348,9 @@ degen_tol : double py::class_> (m, "WorkPlane") .def(py::init(), py::arg("axes")=gp_Ax3(), py::arg("pos")=gp_Ax2d()) - .def_property_readonly("CurrentLocation", &WorkPlane::CurrentLocation) - .def_property_readonly("CurrentDirection", &WorkPlane::CurrentDirection) - .def_property_readonly("StartPnt", &WorkPlane::StartPnt) + .def_property_readonly("cur_loc", &WorkPlane::CurrentLocation) + .def_property_readonly("cur_dir", &WorkPlane::CurrentDirection) + .def_property_readonly("start_pnt", &WorkPlane::StartPnt) .def("MoveTo", &WorkPlane::MoveTo, py::arg("h"), py::arg("v"), "moveto (h,v), and start new wire") .def("Move", &WorkPlane::Move, py::arg("l"), "move 'l' from current position and direction, start new wire") .def("Direction", &WorkPlane::Direction, py::arg("dirh"), py::arg("dirv"), "reset direction to (dirh, dirv)") From 79ff65d224445d415e79f719c495c77f27aa7063 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Fri, 11 Feb 2022 13:56:22 +0100 Subject: [PATCH 12/20] added SIMD-wrapper for 'erf' --- libsrc/core/simd_generic.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/libsrc/core/simd_generic.hpp b/libsrc/core/simd_generic.hpp index 1ad4ea99..b40607a4 100644 --- a/libsrc/core/simd_generic.hpp +++ b/libsrc/core/simd_generic.hpp @@ -604,6 +604,12 @@ namespace ngcore return ngcore::SIMD([a](int i)->double { return log(a[i]); } ); } + using std::erf; + template + NETGEN_INLINE ngcore::SIMD erf (ngcore::SIMD a) { + return ngcore::SIMD([a](int i)->double { return erf(a[i]); } ); + } + using std::pow; template NETGEN_INLINE ngcore::SIMD pow (ngcore::SIMD a, double x) { From 67c031cb7761b7a1f8fcd20d50627f4bf71c17db Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Fri, 11 Feb 2022 15:08:59 +0100 Subject: [PATCH 13/20] fix unused-variable warning --- libsrc/core/profiler.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libsrc/core/profiler.hpp b/libsrc/core/profiler.hpp index c7a60af9..dbe3f1c7 100644 --- a/libsrc/core/profiler.hpp +++ b/libsrc/core/profiler.hpp @@ -147,7 +147,7 @@ namespace ngcore }; }; - + struct TNoTracing{ static constexpr bool do_tracing=false; }; struct TTracing{ static constexpr bool do_tracing=true; }; @@ -163,8 +163,8 @@ namespace ngcore constexpr bool is_timing_type_v = std::is_same_v || std::is_same_v; } - static TNoTracing NoTracing; - static TNoTiming NoTiming; + [[maybe_unused]] static TNoTracing NoTracing; + [[maybe_unused]] static TNoTiming NoTiming; template class Timer From 2ce4412fb2eb1a08325bf773d8f45344f7e52d19 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Fri, 11 Feb 2022 18:19:50 +0100 Subject: [PATCH 14/20] no periodic in bspline surface builder on OCC 7.3 --- libsrc/occ/python_occ_shapes.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 57479a4b..467fb876 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -2204,7 +2204,13 @@ tangents : Dict[int, gp_Vec] points.SetValue(i+1, j+1, gp_Pnt(pnts_unchecked(i, j, 0), pnts_unchecked(i, j, 1), pnts_unchecked(i, j, 2))); GeomAPI_PointsToBSplineSurface builder; +#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 builder.Init(points, approx_type, deg_min, deg_max, continuity, tol, periodic); +#else + if(periodic) + throw Exception("periodic not supported"); + builder.Init(points, approx_type, deg_min, deg_max, continuity, tol); +#endif return BRepBuilderAPI_MakeFace(builder.Surface(), tol).Face(); }, py::arg("points"), @@ -2268,7 +2274,13 @@ degen_tol : double points.SetValue(i+1, j+1, gp_Pnt(pnts_unchecked(i, j, 0), pnts_unchecked(i, j, 1), pnts_unchecked(i, j, 2))); GeomAPI_PointsToBSplineSurface builder; +#if OCC_VERSION_MAJOR>=7 && OCC_VERSION_MINOR>=4 builder.Interpolate(points, approx_type, periodic); +#else + if(periodic) + throw Exception("periodic not supported"); + builder.Interpolate(points, approx_type); +#endif return BRepBuilderAPI_MakeFace(builder.Surface(), degen_tol).Face(); }, py::arg("points"), From 376fe7c694abf322e04c1a2cc6a46c417942a773 Mon Sep 17 00:00:00 2001 From: "mhochsteger@cerbsim.com" Date: Fri, 11 Feb 2022 18:39:48 +0100 Subject: [PATCH 15/20] fix sys import in __main__.py --- python/__main__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/__main__.py b/python/__main__.py index 00ecb5c9..123ae034 100644 --- a/python/__main__.py +++ b/python/__main__.py @@ -1,7 +1,7 @@ -import imp, threading +import imp, threading, sys def handle_arguments(): - import sys, __main__ + import __main__ argv = sys.argv if len(argv)>1 and argv[1].endswith(".py"): with open(argv[1]) as pyfile: From 5d624e3078d7d1e7628830433429de8dcf445c26 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 13 Feb 2022 16:31:55 +0100 Subject: [PATCH 16/20] reading via meshio --- libsrc/meshing/python_mesh.cpp | 83 ++++++++++++++++++++++++++++++++++ libsrc/meshing/topology.cpp | 2 +- python/CMakeLists.txt | 3 +- python/read_meshio.py | 20 ++++++++ 4 files changed, 106 insertions(+), 2 deletions(-) create mode 100644 python/read_meshio.py diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 8881123b..5554e9e8 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -913,6 +913,89 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) { return self.AddFaceDescriptor (fd); }) + + .def ("AddPoints", [](Mesh & self, py::buffer b) + { + py::buffer_info info = b.request(); + if (info.ndim != 2) + throw std::runtime_error("AddPoints needs buffer of dimension 2"); + if (info.format != py::format_descriptor::format()) + throw std::runtime_error("AddPoints needs buffer of type double"); + if (info.strides[0] != sizeof(double)*info.shape[1]) + throw std::runtime_error("AddPoints needs packed array"); + double * ptr = static_cast (info.ptr); + if (info.shape[1]==2) + for (auto i : Range(info.shape[0])) + { + self.AddPoint (Point<3>(ptr[0], ptr[1], 0)); + ptr += 2; + } + if (info.shape[1]==3) + for (auto i : Range(info.shape[0])) + { + self.AddPoint (Point<3>(ptr[0], ptr[1], ptr[2])); + ptr += 3; + } + }) + .def ("AddElements", [](Mesh & self, int dim, int index, py::buffer b, int base) + { + py::buffer_info info = b.request(); + if (info.ndim != 2) + throw std::runtime_error("AddElements needs buffer of dimension 2"); + if (info.format != py::format_descriptor::format()) + throw std::runtime_error("AddPoints needs buffer of type int"); + + int * ptr = static_cast (info.ptr); + if (dim == 2) + { + ELEMENT_TYPE type; + int np = info.shape[1]; + switch (np) + { + case 3: type = TRIG; break; + case 4: type = QUAD; break; + case 6: type = TRIG6; break; + case 8: type = QUAD8; break; + default: + throw Exception("unsupported 2D element with "+ToString(np)+" points"); + } + for (auto i : Range(info.shape[0])) + { + Element2d el(type); + for (int j = 0; j < np;j ++) + el[j] = ptr[j]+PointIndex::BASE-base; + el.SetIndex(index); + self.AddSurfaceElement (el); + ptr += info.strides[0]/sizeof(int); + } + } + if (dim == 3) + { + ELEMENT_TYPE type; + int np = info.shape[1]; + switch (np) + { + case 4: type = TET; break; + /* // have to check ordering of points + case 10: type = TET10; break; + case 8: type = HEX; break; + case 6: type = PRISM; break; + */ + default: + throw Exception("unsupported 3D element with "+ToString(np)+" points"); + } + for (auto i : Range(info.shape[0])) + { + Element el(type); + for (int j = 0; j < np;j ++) + el[j] = ptr[j]+PointIndex::BASE-base; + el.SetIndex(index); + self.AddVolumeElement (el); + ptr += info.strides[0]/sizeof(int); + } + } + + }, py::arg("dim"), py::arg("index"), py::arg("data"), py::arg("base")=0) .def ("DeleteSurfaceElement", [](Mesh & self, SurfaceElementIndex i) diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 3b6e5583..d2a93db2 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -1559,7 +1559,7 @@ namespace netgen } if (cnt_err && ntasks == 1) - cout << cnt_err << " elements are not matching !!!" << endl; + cout << IM(5) << cnt_err << " elements are not matching !!!" << endl; } // NgProfiler::StopTimer (timer2c); diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index c4a648cf..c6c5279f 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -11,7 +11,8 @@ install(FILES ${CMAKE_CURRENT_BINARY_DIR}/config.py ${CMAKE_CURRENT_BINARY_DIR}/version.py __main__.py __init__.py - meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py occ.py read_gmsh.py + meshing.py csg.py geom2d.py stl.py gui.py NgOCC.py occ.py + read_gmsh.py read_meshio.py webgui.py DESTINATION ${NG_INSTALL_DIR_PYTHON}/${NG_INSTALL_SUFFIX} COMPONENT netgen diff --git a/python/read_meshio.py b/python/read_meshio.py new file mode 100644 index 00000000..cbb84bae --- /dev/null +++ b/python/read_meshio.py @@ -0,0 +1,20 @@ +from netgen.meshing import * + +def ReadViaMeshIO(filename): + import meshio + import numpy as np + + # print ("reading via meshio:", filename) + + m = meshio.read(filename) + pts = m.points + + mesh = Mesh(dim=pts.shape[1]) + mesh.AddPoints ( np.asarray(pts, dtype=np.float64) ) # needed for correct little/big endian + + fd = mesh.Add (FaceDescriptor(bc=1)) + for cb in m.cells: + mesh.AddElements(dim=cb.dim, index=1, data=np.asarray(cb.data, dtype=np.int32), base=0) + + return mesh + From ff8708d76bf4bddb9ca45c5f38fc3c27380608af Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 13 Feb 2022 17:29:53 +0100 Subject: [PATCH 17/20] preallocate memory --- libsrc/meshing/python_mesh.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 5554e9e8..30fc1b88 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -916,6 +916,9 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def ("AddPoints", [](Mesh & self, py::buffer b) { + static Timer timer("Mesh::AddPoints"); + RegionTimer reg(timer); + py::buffer_info info = b.request(); if (info.ndim != 2) throw std::runtime_error("AddPoints needs buffer of dimension 2"); @@ -924,6 +927,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) if (info.strides[0] != sizeof(double)*info.shape[1]) throw std::runtime_error("AddPoints needs packed array"); double * ptr = static_cast (info.ptr); + + self.Points().SetAllocSize(self.Points().Size()+info.shape[0]); if (info.shape[1]==2) for (auto i : Range(info.shape[0])) { @@ -939,6 +944,9 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) }) .def ("AddElements", [](Mesh & self, int dim, int index, py::buffer b, int base) { + static Timer timer("Mesh::AddElements"); + RegionTimer reg(timer); + py::buffer_info info = b.request(); if (info.ndim != 2) throw std::runtime_error("AddElements needs buffer of dimension 2"); @@ -959,6 +967,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) default: throw Exception("unsupported 2D element with "+ToString(np)+" points"); } + self.SurfaceElements().SetAllocSize(self.SurfaceElements().Size()+info.shape[0]); for (auto i : Range(info.shape[0])) { Element2d el(type); @@ -984,6 +993,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) default: throw Exception("unsupported 3D element with "+ToString(np)+" points"); } + self.VolumeElements().SetAllocSize(self.VolumeElements().Size()+info.shape[0]); for (auto i : Range(info.shape[0])) { Element el(type); From 6b28275b8862c4b67372573c65acbcda6fcbe3e3 Mon Sep 17 00:00:00 2001 From: Joachim Schoeberl Date: Sun, 13 Feb 2022 19:02:51 +0100 Subject: [PATCH 18/20] double conversion on C++ side --- libsrc/meshing/python_mesh.cpp | 24 ++++++++++++++++++------ python/read_meshio.py | 6 ++++-- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 30fc1b88..38586eb2 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -914,16 +914,23 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) return self.AddFaceDescriptor (fd); }) - .def ("AddPoints", [](Mesh & self, py::buffer b) + .def ("AddPoints", [](Mesh & self, py::buffer b1) { static Timer timer("Mesh::AddPoints"); + static Timer timercast("Mesh::AddPoints - casting"); RegionTimer reg(timer); + + timercast.Start(); + // casting from here: https://github.com/pybind/pybind11/issues/1908 + auto b = b1.cast>(); + timercast.Stop(); py::buffer_info info = b.request(); + // cout << "data format = " << info.format << endl; if (info.ndim != 2) throw std::runtime_error("AddPoints needs buffer of dimension 2"); - if (info.format != py::format_descriptor::format()) - throw std::runtime_error("AddPoints needs buffer of type double"); + // if (info.format != py::format_descriptor::format()) + // throw std::runtime_error("AddPoints needs buffer of type double"); if (info.strides[0] != sizeof(double)*info.shape[1]) throw std::runtime_error("AddPoints needs packed array"); double * ptr = static_cast (info.ptr); @@ -942,16 +949,21 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) ptr += 3; } }) - .def ("AddElements", [](Mesh & self, int dim, int index, py::buffer b, int base) + .def ("AddElements", [](Mesh & self, int dim, int index, py::buffer b1, int base) { static Timer timer("Mesh::AddElements"); + static Timer timercast("Mesh::AddElements casting"); RegionTimer reg(timer); + + timercast.Start(); + auto b = b1.cast>(); + timercast.Stop(); py::buffer_info info = b.request(); if (info.ndim != 2) throw std::runtime_error("AddElements needs buffer of dimension 2"); - if (info.format != py::format_descriptor::format()) - throw std::runtime_error("AddPoints needs buffer of type int"); + // if (info.format != py::format_descriptor::format()) + // throw std::runtime_error("AddPoints needs buffer of type int"); int * ptr = static_cast (info.ptr); if (dim == 2) diff --git a/python/read_meshio.py b/python/read_meshio.py index cbb84bae..5cd8a4f1 100644 --- a/python/read_meshio.py +++ b/python/read_meshio.py @@ -10,11 +10,13 @@ def ReadViaMeshIO(filename): pts = m.points mesh = Mesh(dim=pts.shape[1]) - mesh.AddPoints ( np.asarray(pts, dtype=np.float64) ) # needed for correct little/big endian + # mesh.AddPoints ( np.asarray(pts, dtype=np.float64) ) # needed for correct little/big endian + mesh.AddPoints ( pts ) fd = mesh.Add (FaceDescriptor(bc=1)) for cb in m.cells: - mesh.AddElements(dim=cb.dim, index=1, data=np.asarray(cb.data, dtype=np.int32), base=0) + # mesh.AddElements(dim=cb.dim, index=1, data=np.asarray(cb.data, dtype=np.int32), base=0) + mesh.AddElements(dim=cb.dim, index=1, data=cb.data, base=0) return mesh From 3ee29a1ace5e18947e6ae04f408b1f18bdc5f66b Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Mon, 14 Feb 2022 18:07:00 +0100 Subject: [PATCH 19/20] [occ] overwrite shape property name only if not already set in merge --- libsrc/meshing/basegeom.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 642db838..c381de7c 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -20,7 +20,7 @@ namespace netgen double hpref = 0; // number of hp refinement levels (will be multiplied by factor later) void Merge(const ShapeProperties & prop2) { - if (prop2.name) name = prop2.name; + if (!name && prop2.name) name = prop2.name; if (prop2.col) col = prop2.col; maxh = min2(maxh, prop2.maxh); hpref = max2(hpref, prop2.hpref); From e2040ae95332998ed91ebfb665a57026a462580a Mon Sep 17 00:00:00 2001 From: Christopher Lackner Date: Mon, 14 Feb 2022 18:24:02 +0100 Subject: [PATCH 20/20] [occ] also keep color in merge if already set from shape --- libsrc/meshing/basegeom.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index c381de7c..db82c11c 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -21,7 +21,7 @@ namespace netgen void Merge(const ShapeProperties & prop2) { if (!name && prop2.name) name = prop2.name; - if (prop2.col) col = prop2.col; + if (!col && prop2.col) col = prop2.col; maxh = min2(maxh, prop2.maxh); hpref = max2(hpref, prop2.hpref); }