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 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) { diff --git a/libsrc/meshing/basegeom.hpp b/libsrc/meshing/basegeom.hpp index 642db838..db82c11c 100644 --- a/libsrc/meshing/basegeom.hpp +++ b/libsrc/meshing/basegeom.hpp @@ -20,8 +20,8 @@ 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 (prop2.col) col = prop2.col; + if (!name && prop2.name) name = prop2.name; + if (!col && prop2.col) col = prop2.col; maxh = min2(maxh, prop2.maxh); hpref = max2(hpref, prop2.hpref); } diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 8881123b..38586eb2 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -913,6 +913,111 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) { return self.AddFaceDescriptor (fd); }) + + .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.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])) + { + 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 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"); + + 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"); + } + self.SurfaceElements().SetAllocSize(self.SurfaceElements().Size()+info.shape[0]); + 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"); + } + self.VolumeElements().SetAllocSize(self.VolumeElements().Size()+info.shape[0]); + 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/libsrc/occ/python_occ_shapes.cpp b/libsrc/occ/python_occ_shapes.cpp index 91e94c50..467fb876 100644 --- a/libsrc/occ/python_occ_shapes.cpp +++ b/libsrc/occ/python_occ_shapes.cpp @@ -47,13 +47,19 @@ #include #include #include +#include #include #include #include +#include +#include +#include #include +#include #include -#include #include +#include +#include #include #include #include @@ -230,6 +236,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); @@ -311,6 +331,90 @@ public: return shared_from_this(); } + auto Spline(const std::vector &points, bool periodic, double tol, const std::map &tangents, + bool start_from_localpos) + { + 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()); + + 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]); + } + + 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(); + + const bool closing = periodic || PLast.Distance(startpnt) < 1e-10; + if (startvertex.IsNull()) + startvertex = lastvertex = BRepBuilderAPI_MakeVertex(P13d).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(); + 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 (closing) + Finish(); + + return shared_from_this(); + } + auto ArcTo (double h, double v, const gp_Vec2d t) { gp_Pnt2d P1 = localpos.Location(); @@ -561,8 +665,6 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) .value("SHAPE", TopAbs_SHAPE) .export_values() ; - - py::class_ (m, "TopoDS_Shape") @@ -1597,6 +1699,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); @@ -1719,6 +1837,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 + Maximum 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 { @@ -1879,14 +2091,226 @@ 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", [](const std::vector &points, Approx_ParametrizationType approx_type, int deg_min, + int deg_max, GeomAbs_Shape continuity, double tol) { + 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]); + + GeomAPI_PointsToBSpline builder(hpoints, 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 in 3d. + +Parameters +---------- + +points : List[gp_Pnt] or Tuple[gp_Pnt] + 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(TColgp_HArray1OfPnt) hpoints = new TColgp_HArray1OfPnt(1, points.size()); + for (int i = 0; i < points.size(); i++) + hpoints->SetValue(i+1, points[i]); + + GeomAPI_Interpolate builder(hpoints, periodic, tol); + + 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())) + { + 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 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 in 3d. + +Parameters +---------- + +points : List|Tuple[gp_Pnt] + List (or tuple) of gp_Pnt + +periodic : bool + Whether the result should be periodic + +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"); + + + 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; +#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"), + 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 in the first surface parameter + +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; +#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"), + 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 in the first surface parameter + +degen_tol : double + Tolerance for resolution of degenerate edges + +)delimiter"); m.def("MakeFillet", [](TopoDS_Shape shape, std::vector edges, double r) { @@ -1936,6 +2360,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) py::class_> (m, "WorkPlane") .def(py::init(), py::arg("axes")=gp_Ax3(), py::arg("pos")=gp_Ax2d()) + .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)") @@ -1949,6 +2376,9 @@ DLL_HEADER void ExportNgOCCShapes(py::module &m) 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("start_from_localpos")=true, + "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) { 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/__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: diff --git a/python/read_meshio.py b/python/read_meshio.py new file mode 100644 index 00000000..5cd8a4f1 --- /dev/null +++ b/python/read_meshio.py @@ -0,0 +1,22 @@ +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 + 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=cb.data, base=0) + + return mesh +