diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index df2d341b..683c6db6 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -136,7 +136,7 @@ build_ubuntu_1604: paths: - build/ - src/ - key: "netgen_win32" + key: "netgen_win32_${CI_BUILD_REF_NAME}" test_netgen_win64: <<: *win64 @@ -145,7 +145,7 @@ test_netgen_win64: paths: - build/ - src/ - key: "netgen_win64" + key: "netgen_win64_${CI_BUILD_REF_NAME}" # Linux .template_test_linux: &test_linux diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..fccd5ba0 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "external_dependencies/pybind11"] + path = external_dependencies/pybind11 + url = https://github.com/pybind/pybind11.git diff --git a/CMakeLists.txt b/CMakeLists.txt index edfb7505..13e30bb5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -124,7 +124,7 @@ endmacro() if(WIN32) get_WIN32_WINNT(ver) - add_definitions(-D_WIN32_WINNT=${ver} -DWNT -DWNT_WINDOW) + add_definitions(-D_WIN32_WINNT=${ver} -DWNT -DWNT_WINDOW -DNOMINMAX) set(CMAKE_MFC_FLAG 0) # add_definitions(-DNGINTERFACE_EXPORTS) # add_definitions(-DNGLIB_EXPORTS) @@ -135,7 +135,7 @@ if(WIN32) else(WIN32) # build shared libraries set(NG_LIB_TYPE SHARED) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14") endif(WIN32) if(APPLE) # set(MACOSX_BUNDLE ON) @@ -185,6 +185,17 @@ endif (USE_GUI) ####################################################################### if (USE_PYTHON) + find_path(PYBIND_INCLUDE_DIR pybind11/pybind11.h ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies/pybind11/include) + if( NOT PYBIND_INCLUDE_DIR ) + execute_process(COMMAND git submodule update --init --recursive WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) + find_path(PYBIND_INCLUDE_DIR pybind11/pybind11.h ${CMAKE_CURRENT_SOURCE_DIR}/external_dependencies/pybind11/include) + endif( NOT PYBIND_INCLUDE_DIR ) + if( PYBIND_INCLUDE_DIR ) + message("-- Found Pybind11: ${PYBIND_INCLUDE_DIR}") + else( PYBIND_INCLUDE_DIR ) + message(FATAL_ERROR "Could NOT find pybind11!") + endif( PYBIND_INCLUDE_DIR ) + include_directories(${PYBIND_INCLUDE_DIR}) set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${PROJECT_SOURCE_DIR}/cmake_modules/python") set(PYTHON_VERSION "3" CACHE STRING "") set(Python_ADDITIONAL_VERSIONS 3.5) @@ -198,36 +209,10 @@ if (USE_PYTHON) find_package(PythonInterp ${PYTHON_VERSION} EXACT REQUIRED) find_package(PythonLibs ${PYTHON_VERSION} EXACT REQUIRED) endif( PYTHON_VERSION STREQUAL "3") - set(Boost_FIND_QUIETLY ON) - - foreach( component IN ITEMS python${PYTHON_VERSION_MAJOR} python-py${PYTHON_VERSION_MAJOR}${PYTHON_VERSION_MINOR} python) - string(TOUPPER ${component} component_upper) - if(NOT Boost_LIBRARIES) - find_package(Boost COMPONENTS ${component}) - set(BOOST_PYTHON_LIB ${Boost_${component_upper}_LIBRARY_RELEASE}) - endif(NOT Boost_LIBRARIES) - endforeach() - - if(NOT Boost_LIBRARIES) - message(FATAL_ERROR "Could NOT find Boost-Python" ) - else(NOT Boost_LIBRARIES) - message("-- Found Boost_python: ${BOOST_PYTHON_LIB}") - endif(NOT Boost_LIBRARIES) add_definitions(-DNG_PYTHON) - if(WIN32) - install(FILES ${BOOST_PYTHON_LIB} DESTINATION lib COMPONENT netgen) - install(DIRECTORY ${Boost_INCLUDE_DIR}/boost DESTINATION include COMPONENT netgen) - endif(WIN32) - - if(WIN32 AND ${CMAKE_CXX_COMPILER_ID} STREQUAL "Intel") - # fix problems with boosts config/auto_link.hpp - add_definitions(-DBOOST_LIB_TOOLSET="iw") - endif() - include_directories(${PYTHON_INCLUDE_DIRS}) - include_directories(${Boost_INCLUDE_DIRS}) - set(PYTHON_LIBS "${PYTHON_LIBRARIES};${Boost_LIBRARIES}") + set(PYTHON_LIBS "${PYTHON_LIBRARIES}") execute_process(COMMAND ${PYTHON_EXECUTABLE} -c "from distutils.sysconfig import get_python_lib; print(get_python_lib(1,0,''))" OUTPUT_VARIABLE PYTHON_PACKAGES_INSTALL_DIR OUTPUT_STRIP_TRAILING_WHITESPACE) message("python install path: ${PYTHON_PACKAGES_INSTALL_DIR}") endif (USE_PYTHON) @@ -288,7 +273,7 @@ if(INSTALL_DEPENDENCIES) # find_library(LIBSTDCXX NAMES stdc++.6) # find_library(LIBGCCS NAMES gcc_s.1) # install( FILES "${Boost_LIBRARIES}" ${LIBGOMP} ${LIBSTDCXX} ${LIBGCCS} DESTINATION ${ng_install_dir_lib} COMPONENT netgen ) - install( FILES "${Boost_LIBRARIES}" DESTINATION ${ng_install_dir_lib} COMPONENT netgen ) +# install( FILES "${Boost_LIBRARIES}" DESTINATION ${ng_install_dir_lib} COMPONENT netgen ) get_filename_component(MY_LIB_DIR ${TK_LIBRARY} DIRECTORY) if(APPLE) if(TIX_LIBRARY) @@ -355,8 +340,7 @@ if(UNIX) string(SUBSTRING ${temp} 17 -1 UBUNTU_VERSION) message("ubuntu version: ${UBUNTU_VERSION}") - set(BOOST_PACKAGE "libboost-python${Boost_MAJOR_VERSION}.${Boost_MINOR_VERSION}.${Boost_SUBMINOR_VERSION}") - set(CPACK_DEBIAN_PACKAGE_DEPENDS "python3, ${BOOST_PACKAGE}, libtk8.5, libtcl8.5, tix, libxmu6") + set(CPACK_DEBIAN_PACKAGE_DEPENDS "python3, libtk8.5, libtcl8.5, tix, libxmu6") execute_process(COMMAND dpkg --print-architecture OUTPUT_VARIABLE CPACK_DEBIAN_PACKAGE_ARCHITECTURE OUTPUT_STRIP_TRAILING_WHITESPACE) set(CPACK_DEBIAN_PACKAGE_MAINTAINER "Matthias Hochsteger ") if(USE_MPI) @@ -364,11 +348,7 @@ if(UNIX) set(CPACK_PACKAGE_NAME "${CPACK_PACKAGE_NAME}_mpi") endif(USE_MPI) if(USE_OCC) - if("${UBUNTU_VERSION}" STREQUAL "xenial") - set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf10") - else() - set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf8") - endif() + set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf-dev") endif(USE_OCC) set(CPACK_DEBIAN_PACKAGE_SECTION Science) set(CPACK_DEBIAN_PACKAGE_NAME ${CPACK_PACKAGE_NAME}) diff --git a/external_dependencies/pybind11 b/external_dependencies/pybind11 new file mode 160000 index 00000000..2b92a491 --- /dev/null +++ b/external_dependencies/pybind11 @@ -0,0 +1 @@ +Subproject commit 2b92a49115657712c7dd924248df2286e6143b8d diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 917993a5..fb89cd46 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -1494,7 +1494,7 @@ namespace netgen // virtual VisualScene * GetVisualScene (const NetgenGeometry * geom) const; }; - extern CSGeometry * ParseCSG (istream & istr); + extern CSGeometry * ParseCSG (istream & istr, CSGeometry *instance=nullptr); NetgenGeometry * CSGeometryRegister :: Load (string filename) const { diff --git a/libsrc/csg/csgparser.cpp b/libsrc/csg/csgparser.cpp index c2728bcd..8ac0a6fe 100644 --- a/libsrc/csg/csgparser.cpp +++ b/libsrc/csg/csgparser.cpp @@ -822,11 +822,17 @@ namespace netgen /* Main parsing function for CSG geometry */ - CSGeometry * ParseCSG (istream & istr) + CSGeometry * ParseCSG (istream & istr, CSGeometry * instance=nullptr) { CSGScanner scan(istr); - geom = new CSGeometry; + if (instance) + { + new (instance) CSGeometry; + geom = instance; + } + else + geom = new CSGeometry; scan.ReadNext(); if (scan.GetToken() != TOK_RECO) // keyword 'algebraic3d' diff --git a/libsrc/csg/genmesh.cpp b/libsrc/csg/genmesh.cpp index de45e713..4b8f1713 100644 --- a/libsrc/csg/genmesh.cpp +++ b/libsrc/csg/genmesh.cpp @@ -687,6 +687,8 @@ namespace netgen mparam.grading); mesh -> LoadLocalMeshSize (mparam.meshsizefilename); + for (auto mspnt : mparam.meshsize_points) + mesh -> RestrictLocalH (mspnt.pnt, mspnt.h); } spoints.SetSize(0); diff --git a/libsrc/csg/identify.cpp b/libsrc/csg/identify.cpp index df45773e..ca2d13e6 100644 --- a/libsrc/csg/identify.cpp +++ b/libsrc/csg/identify.cpp @@ -1079,6 +1079,7 @@ void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh) for (int i = 1; i <= mesh.GetNFD(); i++) { + auto & fdi = mesh.GetFaceDescriptor(i); int surfi = mesh.GetFaceDescriptor(i).SurfNr(); if (s1rep != surfi) continue; @@ -1090,12 +1091,21 @@ void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh) for (int j = 1; j <= mesh.GetNFD(); j++) { + auto & fdj = mesh.GetFaceDescriptor(j); int surfj = mesh.GetFaceDescriptor(j).SurfNr(); if (surfi == surfj) continue; if (s2rep != surfj) continue; - + bool have_common = false; + if (fdi.DomainIn() != 0) + if (fdi.DomainIn() == fdj.DomainIn() || fdi.DomainIn() == fdj.DomainOut()) + have_common = true; + if (fdi.DomainOut() != 0) + if (fdi.DomainOut() == fdj.DomainIn() || fdi.DomainOut() == fdj.DomainOut()) + have_common = true; + if (!have_common) continue; + int idok = 1; for (side = 1; side <= 2 && idok; side++) diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 1cc866a8..60db9ee4 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -11,7 +11,6 @@ namespace netgen extern shared_ptr ng_geometry; } -inline void NOOP_Deleter(void *) { ; } // a shadow solid tree using shared pointers. @@ -160,60 +159,103 @@ inline ostream & operator<< (ostream & ost, const SPSolid & sol) namespace netgen { - extern CSGeometry * ParseCSG (istream & istr); + extern CSGeometry * ParseCSG (istream & istr, CSGeometry *instance=nullptr); } -DLL_HEADER void ExportCSG() + +static Transformation<3> global_trafo(Vec<3> (0,0,0)); + +DLL_HEADER void ExportCSG(py::module &m) { - ModuleScope module("csg"); - - bp::class_> ("Point2d", bp::init()) + py::class_(m, "NGDummyArgument") + .def("__bool__", []( NGDummyArgument &self ) { return false; } ) + ; + + py::class_> (m, "Point2d") + .def(py::init()) .def ("__str__", &ToString>) - .def(bp::self-bp::self) - .def(bp::self+Vec<2>()) - .def(bp::self-Vec<2>()) + .def(py::self-py::self) + .def(py::self+Vec<2>()) + .def(py::self-Vec<2>()) ; - bp::class_> ("Point3d", bp::init()) + py::class_> (m, "Point3d") + .def(py::init()) .def ("__str__", &ToString>) - .def(bp::self-bp::self) - .def(bp::self+Vec<3>()) - .def(bp::self-Vec<3>()) + .def(py::self-py::self) + .def(py::self+Vec<3>()) + .def(py::self-Vec<3>()) ; - bp::def ("Pnt", FunctionPointer - ([](double x, double y, double z) { return Point<3>(x,y,z); })); - bp::def ("Pnt", FunctionPointer + m.def ("Pnt", FunctionPointer + ([](double x, double y, double z) { return global_trafo(Point<3>(x,y,z)); })); + m.def ("Pnt", FunctionPointer ([](double x, double y) { return Point<2>(x,y); })); - bp::class_> ("Vec2d", bp::init()) + + m.def ("Pnt", FunctionPointer + ([](double x, double y, double z) { return Point<3>(x,y,z); })); + m.def ("Pnt", FunctionPointer + ([](double x, double y) { return Point<2>(x,y); })); + + m.def ("SetTransformation", FunctionPointer + ([](int dir, double angle) + { + if (dir > 0) + global_trafo.SetAxisRotation (dir, angle*M_PI/180); + else + global_trafo = Transformation<3> (Vec<3>(0,0,0)); + }), + py::arg("dir")=int(0), py::arg("angle")=int(0)); + + py::class_> (m, "Vec2d") + .def(py::init()) .def ("__str__", &ToString>) - .def(bp::self+bp::self) - .def(bp::self-bp::self) - .def(-bp::self) - .def(double()*bp::self) + .def(py::self+py::self) + .def(py::self-py::self) + .def(-py::self) + .def(double()*py::self) .def("Norm", &Vec<2>::Length) ; - bp::class_> ("Vec3d", bp::init()) + py::class_> (m, "Vec3d") + .def(py::init()) .def ("__str__", &ToString>) - .def(bp::self+bp::self) - .def(bp::self-bp::self) - .def(-bp::self) - .def(double()*bp::self) + .def(py::self+py::self) + .def(py::self-py::self) + .def(-py::self) + .def(double()*py::self) .def("Norm", &Vec<3>::Length) ; - bp::def ("Vec", FunctionPointer - ([] (double x, double y, double z) { return Vec<3>(x,y,z); })); - bp::def ("Vec", FunctionPointer + m.def ("Vec", FunctionPointer + ([] (double x, double y, double z) { return global_trafo(Vec<3>(x,y,z)); })); + m.def ("Vec", FunctionPointer ([] (double x, double y) { return Vec<2>(x,y); })); -#if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100) - bp::register_ptr_to_python>(); -#endif - bp::class_, boost::noncopyable> ("Solid", bp::no_init) + py::class_> (m, "SplineCurve2d") + .def(py::init<>()) + .def ("AddPoint", FunctionPointer + ([] (SplineGeometry<2> & self, double x, double y) + { + self.geompoints.Append (GeomPoint<2> (Point<2> (x,y))); + return self.geompoints.Size()-1; + })) + .def ("AddSegment", FunctionPointer + ([] (SplineGeometry<2> & self, int i1, int i2) + { + self.splines.Append (new LineSeg<2> (self.geompoints[i1], self.geompoints[i2])); + })) + .def ("AddSegment", FunctionPointer + ([] (SplineGeometry<2> & self, int i1, int i2, int i3) + { + self.splines.Append (new SplineSeg3<2> (self.geompoints[i1], self.geompoints[i2], self.geompoints[i3])); + })) + ; + + + py::class_> (m, "Solid") .def ("__str__", &ToString) .def ("__add__", FunctionPointer( [] ( shared_ptr self, shared_ptr other) { return make_shared (SPSolid::UNION, self, other); })) .def ("__mul__", FunctionPointer( [] ( shared_ptr self, shared_ptr other) { return make_shared (SPSolid::SECTION, self, other); })) @@ -231,100 +273,105 @@ DLL_HEADER void ExportCSG() .def ("mat", FunctionPointer([](shared_ptr & self, string mat) -> shared_ptr { self->SetMaterial(mat); return self; })) .def ("mat", &SPSolid::GetMaterial) - .def("col", FunctionPointer([](shared_ptr & self, bp::list rgb) -> shared_ptr + .def("col", FunctionPointer([](shared_ptr & self, py::list rgb) -> shared_ptr { - bp::extract red(rgb[0]); - bp::extract green(rgb[1]); - bp::extract blue(rgb[2]); + py::extract red(rgb[0]); + py::extract green(rgb[1]); + py::extract blue(rgb[2]); self->SetColor(red(),green(),blue()); return self; })) .def("transp", FunctionPointer([](shared_ptr & self)->shared_ptr < SPSolid > { self->SetTransparent(); return self; })) ; - bp::def ("Sphere", FunctionPointer([](Point<3> c, double r) + m.def ("Sphere", FunctionPointer([](Point<3> c, double r) { Sphere * sp = new Sphere (c, r); Solid * sol = new Solid (sp); return make_shared (sol); })); - bp::def ("Plane", FunctionPointer([](Point<3> p, Vec<3> n) + m.def ("Plane", FunctionPointer([](Point<3> p, Vec<3> n) { Plane * sp = new Plane (p,n); Solid * sol = new Solid (sp); return make_shared (sol); })); - bp::def ("Cone", FunctionPointer([](Point<3> a, Point<3> b, double ra, double rb) + m.def ("Cone", FunctionPointer([](Point<3> a, Point<3> b, double ra, double rb) { Cone * cyl = new Cone (a, b, ra, rb); Solid * sol = new Solid (cyl); return make_shared (sol); })); - bp::def ("Cylinder", FunctionPointer([](Point<3> a, Point<3> b, double r) + m.def ("Cylinder", FunctionPointer([](Point<3> a, Point<3> b, double r) { Cylinder * cyl = new Cylinder (a, b, r); Solid * sol = new Solid (cyl); return make_shared (sol); })); - bp::def ("OrthoBrick", FunctionPointer([](Point<3> p1, Point<3> p2) + m.def ("OrthoBrick", FunctionPointer([](Point<3> p1, Point<3> p2) { OrthoBrick * brick = new OrthoBrick (p1,p2); Solid * sol = new Solid (brick); return make_shared (sol); })); + m.def ("Revolution", FunctionPointer([](Point<3> p1, Point<3> p2, + const SplineGeometry<2> & spline) + { + Revolution * rev = new Revolution (p1, p2, spline); + Solid * sol = new Solid(rev); + return make_shared (sol); + })); - bp::def ("Or", FunctionPointer([](shared_ptr s1, shared_ptr s2) + m.def ("Or", FunctionPointer([](shared_ptr s1, shared_ptr s2) { return make_shared (SPSolid::UNION, s1, s2); })); - bp::def ("And", FunctionPointer([](shared_ptr s1, shared_ptr s2) + m.def ("And", FunctionPointer([](shared_ptr s1, shared_ptr s2) { return make_shared (SPSolid::SECTION, s1, s2); })); - bp::class_, boost::noncopyable> ("CSGeometry") - .def("__init__", bp::make_constructor (FunctionPointer - ([](const string & filename) + py::class_> (m, "CSGeometry") + .def(py::init<>()) + .def("__init__", + [](CSGeometry *instance, const string & filename) { cout << "load geometry"; ifstream ist(filename); - shared_ptr geom(ParseCSG(ist)); - geom -> FindIdenticSurfaces(1e-8 * geom->MaxSize()); - return geom; - }))) - - .def("__init__", bp::make_constructor (FunctionPointer - ([](const bp::list & solidlist) + ParseCSG(ist, instance); + instance -> FindIdenticSurfaces(1e-8 * instance->MaxSize()); + }) + .def("__init__", + [](CSGeometry *instance, const py::list & solidlist) { cout << "csg from list"; - auto geom = make_shared(); + new (instance) CSGeometry(); for (int i = 0; i < len(solidlist); i++) { - bp::object obj = solidlist[i]; + py::object obj = solidlist[i]; cout << "obj " << i << endl; - bp::extract> solid(solidlist[i]); + py::extract> solid(solidlist[i]); if(solid.check()) { cout << "its a solid" << endl; - solid()->AddSurfaces (*geom); + solid()->AddSurfaces (*instance); solid()->GiveUpOwner(); - int tlonr = geom->SetTopLevelObject (solid()->GetSolid()); - geom->GetTopLevelObject(tlonr) -> SetMaterial(solid()->GetMaterial()); + int tlonr = instance->SetTopLevelObject (solid()->GetSolid()); + instance->GetTopLevelObject(tlonr) -> SetMaterial(solid()->GetMaterial()); } } - geom -> FindIdenticSurfaces(1e-8 * geom->MaxSize()); - return geom; - }))) + instance -> FindIdenticSurfaces(1e-8 * instance->MaxSize()); + }) .def("Save", FunctionPointer([] (CSGeometry & self, string filename) { cout << "save geometry to file " << filename << endl; self.Save (filename); })) - .def("Add", FunctionPointer - ([] (CSGeometry & self, shared_ptr solid, bp::list bcmod) + .def("Add", + [] (CSGeometry & self, shared_ptr solid, py::list bcmod) { solid->AddSurfaces (self); solid->GiveUpOwner(); @@ -334,15 +381,15 @@ DLL_HEADER void ExportCSG() self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent()); // bcmod is list of tuples ( solid, bcnr ) - for (int i = 0; i < bp::len(bcmod); i++) + for (int i = 0; i < py::len(bcmod); i++) { - bp::tuple tup = bp::extract (bcmod[i]) (); - auto mod_solid = bp::extract> (tup[0]) (); + py::tuple tup = py::extract (bcmod[i]) (); + auto mod_solid = py::extract> (tup[0]) (); int mod_nr = -1; string * bcname = nullptr; - bp::object val = tup[1]; - if (bp::extract(val).check()) mod_nr = bp::extract (val)(); - if (bp::extract(val).check()) bcname = new string ( bp::extract (val)()); + py::object val = tup[1]; + if (py::extract(val).check()) mod_nr = py::extract (val)(); + if (py::extract(val).check()) bcname = new string ( py::extract (val)()); Array si; mod_solid -> GetSolid() -> GetSurfaceIndices (si); @@ -360,9 +407,8 @@ DLL_HEADER void ExportCSG() delete bcname; } return tlonr; - - }), - (bp::arg("self"), bp::arg("solid"), bp::arg("bcmod")=bp::list()) + }, + py::arg("solid"), py::arg("bcmod")=py::list() ) .def("AddSurface", FunctionPointer @@ -379,13 +425,13 @@ DLL_HEADER void ExportCSG() self.GetTopLevelObject(tlonr) -> SetRGB(solid->GetRed(),solid->GetGreen(),solid->GetBlue()); self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent()); }), - (bp::arg("self"), bp::arg("surface"), bp::arg("solid")) + py::arg("surface"), py::arg("solid") ) .def("CloseSurfaces", FunctionPointer - ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, bp::list aslices ) + ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, py::list aslices ) { Array si1, si2; s1->GetSolid()->GetSurfaceIndices (si1); @@ -397,15 +443,15 @@ DLL_HEADER void ExportCSG() try { - int n = bp::len(aslices); + int n = py::len(aslices); Array slices(n); for(int i=0; i(aslices[i])(); + slices[i]= py::extract(aslices[i])(); } flags.SetFlag("slices", slices); } - catch( bp::error_already_set const & ) { + catch( py::error_already_set const & ) { cout << "caught python error:" << endl; PyErr_Print(); } @@ -418,7 +464,7 @@ DLL_HEADER void ExportCSG() domain, flags)); }), - (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2"), bp::arg("slices")) + py::arg("solid1"), py::arg("solid2"), py::arg("slices") ) .def("CloseSurfaces", FunctionPointer ([] (CSGeometry & self, shared_ptr s1, shared_ptr s2, int reflevels) @@ -438,7 +484,7 @@ DLL_HEADER void ExportCSG() domain, flags)); }), - (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2"), bp::arg("reflevels")=2) + py::arg("solid1"), py::arg("solid2"), py::arg("reflevels")=2 ) .def("PeriodicSurfaces", FunctionPointer @@ -453,7 +499,7 @@ DLL_HEADER void ExportCSG() (self.GetNIdentifications()+1, self, self.GetSurface (si1[0]), self.GetSurface (si2[0]))); }), - (bp::arg("self"), bp::arg("solid1"), bp::arg("solid2")) + py::arg("solid1"), py::arg("solid2") ) .def("GetTransparent", FunctionPointer @@ -461,14 +507,14 @@ DLL_HEADER void ExportCSG() { return self.GetTopLevelObject(tlonr)->GetTransparent(); }), - (bp::arg("self"), bp::arg("tlonr")) + py::arg("tlonr") ) .def("SetTransparent", FunctionPointer ([] (CSGeometry & self, int tlonr, bool transparent) { self.GetTopLevelObject(tlonr)->SetTransparent(transparent); }), - (bp::arg("self"), bp::arg("tlonr"), bp::arg("transparent")) + py::arg("tlonr"), py::arg("transparent") ) .def("GetVisible", FunctionPointer @@ -476,26 +522,34 @@ DLL_HEADER void ExportCSG() { return self.GetTopLevelObject(tlonr)->GetVisible(); }), - (bp::arg("self"), bp::arg("tlonr")) + py::arg("tlonr") ) .def("SetVisible", FunctionPointer ([] (CSGeometry & self, int tlonr, bool visible) { self.GetTopLevelObject(tlonr)->SetVisible(visible); }), - (bp::arg("self"), bp::arg("tlonr"), bp::arg("visible")) + py::arg("tlonr"), py::arg("visible") ) .def("SetBoundingBox", FunctionPointer ([] (CSGeometry & self, Point<3> pmin, Point<3> pmax) { self.SetBoundingBox(Box<3> (pmin, pmax)); }), - (bp::arg("self"), bp::arg("pmin"), bp::arg("pmax")) + py::arg("pmin"), py::arg("pmax") ) - .add_property ("ntlo", &CSGeometry::GetNTopLevelObjects) + .def("Draw", FunctionPointer + ([] (shared_ptr self) + { + self->FindIdenticSurfaces(1e-6); + self->CalcTriangleApproximation(0.01, 20); + ng_geometry = self; + }) + ) + .def_property_readonly ("ntlo", &CSGeometry::GetNTopLevelObjects) ; - bp::def("GenerateMesh", FunctionPointer + m.def("GenerateMesh", FunctionPointer ([](shared_ptr geo, MeshingParameters & param) { auto dummy = make_shared(); @@ -515,7 +569,7 @@ DLL_HEADER void ExportCSG() })) ; - bp::def("Save", FunctionPointer + m.def("Save", FunctionPointer ([](const Mesh & self, const string & filename, const CSGeometry & geom) { ostream * outfile; @@ -533,7 +587,7 @@ DLL_HEADER void ExportCSG() - bp::def("ZRefinement", FunctionPointer + m.def("ZRefinement", FunctionPointer ([](Mesh & mesh, CSGeometry & geom) { ZRefinementOptions opt; @@ -543,14 +597,10 @@ DLL_HEADER void ExportCSG() ; } - - - - -BOOST_PYTHON_MODULE(libcsg) { - ExportCSG(); +PYBIND11_PLUGIN(libcsg) { + py::module m("csg", "pybind csg"); + ExportCSG(m); + return m.ptr(); } - - #endif diff --git a/libsrc/csg/revolution.cpp b/libsrc/csg/revolution.cpp index 1b48676e..85b49bd8 100644 --- a/libsrc/csg/revolution.cpp +++ b/libsrc/csg/revolution.cpp @@ -156,25 +156,42 @@ namespace netgen { if(spline_coefficient.Size() == 0) spline->GetCoeff(spline_coefficient); + if(spline_coefficient_shifted.Size() == 0) + spline->GetCoeff(spline_coefficient_shifted, spline->StartPI()); Point<2> p; CalcProj(point,p); - return spline_coefficient(0)*p(0)*p(0) + spline_coefficient(1)*p(1)*p(1) + double val = spline_coefficient(0)*p(0)*p(0) + spline_coefficient(1)*p(1)*p(1) + spline_coefficient(2)*p(0)*p(1) + spline_coefficient(3)*p(0) + spline_coefficient(4)*p(1) + spline_coefficient(5); + + Vec<2> pr = p-spline->StartPI(); + + + // cout << "spline_coefficinet = " << spline_coefficient << endl; + // cout << "shifted = " << spline_coefficient_shifted << endl; + + double val2 = spline_coefficient_shifted(0)*pr(0)*pr(0) + spline_coefficient_shifted(1)*pr(1)*pr(1) + + spline_coefficient_shifted(2)*pr(0)*pr(1) + spline_coefficient_shifted(3)*pr(0) + + spline_coefficient_shifted(4)*pr(1) + spline_coefficient_shifted(5); + + // cout << "val = " << val << " =?= " << val2 << endl; + return val2; } void RevolutionFace :: CalcGradient (const Point<3> & point, Vec<3> & grad) const { if(spline_coefficient.Size() == 0) spline->GetCoeff(spline_coefficient); + if(spline_coefficient_shifted.Size() == 0) + spline->GetCoeff(spline_coefficient_shifted, spline->StartPI()); Vec<3> point_minus_p0 = point-p0; Point<2> p; CalcProj0(point_minus_p0,p); - + /* const double dFdxbar = 2.*spline_coefficient(0)*p(0) + spline_coefficient(2)*p(1) + spline_coefficient(3); if(fabs(p(1)) > 1e-10) @@ -193,6 +210,27 @@ namespace netgen grad(2) = dFdxbar*v_axis(2); //(*testout) << "grad2("< pr = p-spline->StartPI(); + const double dFdxbar = 2.*spline_coefficient_shifted(0)*pr(0) + spline_coefficient_shifted(2)*pr(1) + spline_coefficient_shifted(3); + + if(fabs(p(1)) > 1e-10) + { + const double dFdybar = 2.*spline_coefficient_shifted(1)*pr(1) + spline_coefficient_shifted(2)*pr(0) + spline_coefficient_shifted(4); + + grad(0) = dFdxbar*v_axis(0) + dFdybar * ( point_minus_p0(0)-v_axis(0)*p(0) )/p(1); + grad(1) = dFdxbar*v_axis(1) + dFdybar * ( point_minus_p0(1)-v_axis(1)*p(0) )/p(1); + grad(2) = dFdxbar*v_axis(2) + dFdybar * ( point_minus_p0(2)-v_axis(2)*p(0) )/p(1); + //(*testout) << "grad1("<> - ("VisualSceneGeometry", bp::no_init) + py::class_> + (m, "VisualSceneGeometry") .def("Draw", &VisualSceneGeometry::DrawScene) ; - bp::def("SetBackGroundColor", &VisualSceneGeometry::SetBackGroundColor); + m.def("SetBackGroundColor", &VisualSceneGeometry::SetBackGroundColor); - bp::def("VS", FunctionPointer - ([](CSGeometry & geom) + m.def("VS", + [](CSGeometry & geom) { geom.FindIdenticSurfaces(1e-6); geom.CalcTriangleApproximation(0.01, 20); @@ -514,17 +510,19 @@ DLL_HEADER void ExportCSGVis() vs->SetGeometry(&geom); return vs; - })); + }); - bp::def("MouseMove", FunctionPointer - ([](VisualSceneGeometry &vsgeom, int oldx, int oldy, int newx, int newy, char mode) + m.def("MouseMove", + [](VisualSceneGeometry &vsgeom, int oldx, int oldy, int newx, int newy, char mode) { vsgeom.MouseMove(oldx, oldy, newx, newy, mode); - })); + }); } -BOOST_PYTHON_MODULE(libcsgvis) -{ - ExportCSGVis(); + +PYBIND11_PLUGIN(libcsgvis) { + py::module m("csg", "pybind csg"); + ExportCSGVis(m); + return m.ptr(); } #endif diff --git a/libsrc/general/hashtabl.hpp b/libsrc/general/hashtabl.hpp index 18e19fce..b856bbc6 100644 --- a/libsrc/general/hashtabl.hpp +++ b/libsrc/general/hashtabl.hpp @@ -779,7 +779,7 @@ protected: mask = size-1; // cout << "mask = " << mask << endl; invalid = -1; - for (int i = 0; i < size; i++) + for (size_t i = 0; i < size; i++) hash[i].I1() = invalid; } diff --git a/libsrc/general/mpi_interface.cpp b/libsrc/general/mpi_interface.cpp index 57d69ab9..36df5542 100644 --- a/libsrc/general/mpi_interface.cpp +++ b/libsrc/general/mpi_interface.cpp @@ -16,17 +16,17 @@ namespace netgen void MyMPI_SendCmd (const char * cmd) { - char buf[100]; + char buf[10000]; strcpy (buf, cmd); // MPI_Bcast (&buf, 100, MPI_CHAR, 0, MPI_COMM_WORLD); for (int dest = 1; dest < ntasks; dest++) - MPI_Bsend( &buf, 100, MPI_CHAR, dest, MPI_TAG_CMD, MPI_COMM_WORLD); + MPI_Send( &buf, 10000, MPI_CHAR, dest, MPI_TAG_CMD, MPI_COMM_WORLD); } string MyMPI_RecvCmd () { - char buf[100]; + char buf[10000]; // MPI_Bcast (&buf, 100, MPI_CHAR, 0, MPI_COMM_WORLD); // VT_OFF(); @@ -44,7 +44,7 @@ namespace netgen while (!flag); // VT_ON(); - MPI_Recv( &buf, 100, MPI_CHAR, 0, MPI_TAG_CMD, MPI_COMM_WORLD, &status); + MPI_Recv( &buf, 10000, MPI_CHAR, 0, MPI_TAG_CMD, MPI_COMM_WORLD, &status); return string(buf); } diff --git a/libsrc/general/mysimd.hpp b/libsrc/general/mysimd.hpp index 621ea5a3..bc00adae 100644 --- a/libsrc/general/mysimd.hpp +++ b/libsrc/general/mysimd.hpp @@ -54,8 +54,19 @@ namespace netgen #ifdef __AVX__ + template + class AlignedAlloc + { + public: + void * operator new (size_t s, void *p) { return p; } + void * operator new (size_t s) { return _mm_malloc(s, alignof(T)); } + void * operator new[] (size_t s) { return _mm_malloc(s, alignof(T)); } + void operator delete (void * p) { _mm_free(p); } + void operator delete[] (void * p) { _mm_free(p); } + }; + template<> - class alignas(32) SIMD + class alignas(32) SIMD : public AlignedAlloc> { __m256d data; @@ -158,6 +169,10 @@ namespace netgen #else + // it's only a dummy without AVX + template + class AlignedAlloc { ; }; + template<> class SIMD { diff --git a/libsrc/general/ngpython.hpp b/libsrc/general/ngpython.hpp index 2ccad390..7c2d98a5 100644 --- a/libsrc/general/ngpython.hpp +++ b/libsrc/general/ngpython.hpp @@ -1,36 +1,42 @@ #ifdef NG_PYTHON -#include -namespace bp = boost::python; +#include +#include +namespace py = pybind11; #include +#include + +namespace pybind11 { +template +bool CheckCast( py::handle obj ) { + try{ + obj.cast(); + return true; + } + catch (py::cast_error &e) { + return false; + } +} + + +template +struct extract +{ + py::handle obj; + extract( py::handle aobj ) : obj(aobj) {} + + bool check() { return CheckCast(obj); } + T operator()() { return obj.cast(); } +}; +} + +struct NGDummyArgument {}; + +inline void NOOP_Deleter(void *) { ; } namespace netgen { - class ModuleScope { - bp::scope *local_scope; - public: - ModuleScope(const std::string name) : local_scope(nullptr) { - std::string nested_name = name; - if (bp::scope()) - nested_name = bp::extract(bp::scope().attr("__name__") + "." + name); - - bp::object module(bp::handle<>(bp::borrowed(PyImport_AddModule(nested_name.c_str())))); - - std::cout << "exporting " << nested_name << std::endl; - bp::object parent = bp::scope() ? bp::scope() : bp::import("__main__"); - parent.attr(name.c_str()) = module; - - local_scope = new bp::scope(module); - } - - ~ModuleScope() { - if (local_scope) - delete (local_scope); - } - - }; - ////////////////////////////////////////////////////////////////////// // Lambda to function pointer conversion template diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index 69378a3e..29bed482 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -109,7 +109,12 @@ namespace netgen double curv = fabs(first(0)*second(1)-first(1)*second(0)) / pow(first.Length(), 3); return curv; } - + + virtual bool InConvexHull (Point<2> p, double eps) const + { + return seg.InConvexHull (p, eps); + } + }; diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index 8c3f1f42..7dfc1b4c 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -1,13 +1,11 @@ #ifdef NG_PYTHON -#include #include <../general/ngpython.hpp> #include #include using namespace netgen; -namespace bp = boost::python; namespace netgen { @@ -15,24 +13,22 @@ namespace netgen } -DLL_HEADER void ExportGeom2d() +DLL_HEADER void ExportGeom2d(py::module &m) { - ModuleScope module("geom2d"); - bp::class_, boost::noncopyable> - ("SplineGeometry", + py::class_> + (m, "SplineGeometry", "a 2d boundary representation geometry model by lines and splines") - .def("__init__", bp::make_constructor - (FunctionPointer - ([](const string & filename) + .def(py::init<>()) + .def("__init__", + [](SplineGeometry2d *instance, const string & filename) { cout << "load geometry"; ifstream ist(filename); - auto geom = make_shared(); - geom->Load (filename.c_str()); - ng_geometry = geom; - return geom; - }))) + new (instance) SplineGeometry2d(); + instance->Load (filename.c_str()); + ng_geometry = shared_ptr(instance, NOOP_Deleter); + }) .def("Load",&SplineGeometry2d::Load) .def("AppendPoint", FunctionPointer @@ -47,16 +43,17 @@ DLL_HEADER void ExportGeom2d() self.geompoints.Append(gp); return self.geompoints.Size()-1; }), - (bp::arg("self"), bp::arg("x"), bp::arg("y"), bp::arg("maxh") = 1e99, bp::arg("hpref")=false)) - .def("Append", FunctionPointer([](SplineGeometry2d &self, bp::list segment, int leftdomain, int rightdomain, bp::object bc, double maxh, bool hpref) + py::arg("x"), py::arg("y"), py::arg("maxh") = 1e99, py::arg("hpref")=false) + .def("Append", FunctionPointer([](SplineGeometry2d &self, py::list segment, int leftdomain, int rightdomain, + py::object bc, py::object copy, double maxh, bool hpref) { - bp::extract segtype(segment[0]); + py::extract segtype(segment[0]); SplineSegExt * seg; if (segtype().compare("line") == 0) { - bp::extract point_index1(segment[1]); - bp::extract point_index2(segment[2]); + py::extract point_index1(segment[1]); + py::extract point_index2(segment[2]); //point_index1.check() LineSeg<2> * l = new LineSeg<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2())); @@ -64,9 +61,9 @@ DLL_HEADER void ExportGeom2d() } else if (segtype().compare("spline3") == 0) { - bp::extract point_index1(segment[1]); - bp::extract point_index2(segment[2]); - bp::extract point_index3(segment[3]); + py::extract point_index1(segment[1]); + py::extract point_index2(segment[2]); + py::extract point_index3(segment[3]); SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2()), self.GetPoint(point_index3())); seg = new SplineSegExt(*seg3); @@ -82,11 +79,14 @@ DLL_HEADER void ExportGeom2d() seg->hpref_right = hpref; seg->reffak = 1; seg->copyfrom = -1; - if (bp::extract(bc).check()) - seg->bc = bp::extract(bc)(); - else if (bp::extract(bc).check()) + if (py::extract(copy).check()) + seg->copyfrom = py::extract(copy)()+1; + + if (py::extract(bc).check()) + seg->bc = py::extract(bc)(); + else if (py::extract(bc).check()) { - string bcname = bp::extract(bc)(); + string bcname = py::extract(bc)(); int bcnum = self.GetBCNumber(bcname); if (bcnum == 0) bcnum = self.AddBCName(bcname); @@ -95,25 +95,26 @@ DLL_HEADER void ExportGeom2d() else seg->bc = self.GetNSplines()+1; self.AppendSegment(seg); - }), (bp::arg("self"), bp::arg("point_indices"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0, - bp::arg("bc")=bp::object(), bp::arg("maxh")=1e99, bp::arg("hpref")=false - )) + return self.GetNSplines()-1; + }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = py::int_(0), + py::arg("bc")=NGDummyArgument(), py::arg("copy")=NGDummyArgument(), py::arg("maxh")=1e99, py::arg("hpref")=false + ) - .def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, bp::list point_indices, int leftdomain, int rightdomain) + .def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, py::list point_indices, int leftdomain, int rightdomain) { - int npts = bp::len(point_indices); + int npts = py::len(point_indices); SplineSegExt * seg; - //int a = bp::extract(point_indices[0]); + //int a = py::extract(point_indices[0]); if (npts == 2) { - LineSeg<2> * l = new LineSeg<2>(self.GetPoint(bp::extract(point_indices[0])), self.GetPoint(bp::extract(point_indices[1]))); + LineSeg<2> * l = new LineSeg<2>(self.GetPoint(py::extract(point_indices[0])()), self.GetPoint(py::extract(point_indices[1])())); seg = new SplineSegExt(*l); } else if (npts == 3) { - SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(bp::extract(point_indices[0])), self.GetPoint(bp::extract(point_indices[1])), self.GetPoint(bp::extract(point_indices[2]))); + SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(py::extract(point_indices[0])()), self.GetPoint(py::extract(point_indices[1])()), self.GetPoint(py::extract(point_indices[2])())); seg = new SplineSegExt(*seg3); } @@ -123,7 +124,7 @@ DLL_HEADER void ExportGeom2d() seg->reffak = 1; seg->copyfrom = -1; self.AppendSegment(seg); - }), (bp::arg("self"), bp::arg("point_indices"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0) ) + }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = py::int_(0)) //.def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, int point_index1, int point_index2)//, int leftdomain, int rightdomain) // { // LineSeg<2> * l = new LineSeg<2>(self.GetPoint(point_index1), self.GetPoint(point_index2)); @@ -135,7 +136,7 @@ DLL_HEADER void ExportGeom2d() // seg->copyfrom = -1; // self.AppendSegment(seg); - // }))//, (bp::arg("self"), bp::arg("point_index1"), bp::arg("point_index2"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0) ) + // }))//, (py::arg("self"), py::arg("point_index1"), py::arg("point_index2"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0) ) //.def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, int point_index1, int point_index2, int point_index3)//, int leftdomain, int rightdomain) // { // SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1), self.GetPoint(point_index2), self.GetPoint(point_index3)); @@ -146,7 +147,7 @@ DLL_HEADER void ExportGeom2d() // seg->reffak = 1; // seg->copyfrom = -1; // self.AppendSegment(seg); - // }))//, (bp::arg("self"), bp::arg("point_index1"), bp::arg("point_index2"), bp::arg("point_index3"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0 ) ) + // }))//, (py::arg("self"), py::arg("point_index1"), py::arg("point_index2"), py::arg("point_index3"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0 ) ) .def("SetMaterial", &SplineGeometry2d::SetMaterial) @@ -161,22 +162,22 @@ DLL_HEADER void ExportGeom2d() Box<2> box(self.GetBoundingBox()); double xdist = box.PMax()(0) - box.PMin()(0); double ydist = box.PMax()(1) - box.PMin()(1); - bp::tuple xlim = bp::make_tuple(box.PMin()(0) - 0.1*xdist, box.PMax()(0) + 0.1*xdist); - bp::tuple ylim = bp::make_tuple(box.PMin()(1) - 0.1*ydist, box.PMax()(1) + 0.1*ydist); + py::tuple xlim = py::make_tuple(box.PMin()(0) - 0.1*xdist, box.PMax()(0) + 0.1*xdist); + py::tuple ylim = py::make_tuple(box.PMin()(1) - 0.1*ydist, box.PMax()(1) + 0.1*ydist); - bp::list xpoints, ypoints; + py::list xpoints, ypoints; for (int i = 0; i < self.splines.Size(); i++) { - bp::list xp, yp; + py::list xp, yp; if (self.splines[i]->GetType().compare("line")==0) { GeomPoint<2> p1 = self.splines[i]->StartPI(); GeomPoint<2> p2 = self.splines[i]->EndPI(); - xp.append(p1(0)); - xp.append(p2(0)); - yp.append(p1(1)); - yp.append(p2(1)); + xp.append(py::cast(p1(0))); + xp.append(py::cast(p2(0))); + yp.append(py::cast(p1(1))); + yp.append(py::cast(p2(1))); } else if (self.splines[i]->GetType().compare("spline3")==0) { @@ -186,37 +187,37 @@ DLL_HEADER void ExportGeom2d() for (int j = 0; j <= n; j++) { GeomPoint<2> point = self.splines[i]->GetPoint(j*1./n); - xp.append(point(0)); - yp.append(point(1)); + xp.append(py::cast(point(0))); + yp.append(py::cast(point(1))); } } else { cout << "spline is neither line nor spline3" << endl; } - xpoints.append(xp); - ypoints.append(yp); + xpoints.append(py::cast(xp)); + ypoints.append(py::cast(yp)); } - return bp::tuple(bp::make_tuple(xlim, ylim, xpoints, ypoints)); + return py::tuple(py::make_tuple(xlim, ylim, xpoints, ypoints)); })) .def("PointData", FunctionPointer([](SplineGeometry2d &self) { - bp::list xpoints, ypoints, pointindex; + py::list xpoints, ypoints, pointindex; for (int i = 0; i < self.geompoints.Size(); i++) { - pointindex.append(i); - xpoints.append(self.geompoints[i][0]); - ypoints.append(self.geompoints[i][1]); + pointindex.append(py::cast(i)); + xpoints.append(py::cast(self.geompoints[i][0])); + ypoints.append(py::cast(self.geompoints[i][1])); } - return bp::tuple(bp::make_tuple(xpoints, ypoints, pointindex)); + return py::tuple(py::make_tuple(xpoints, ypoints, pointindex)); })) .def("SegmentData", FunctionPointer([](SplineGeometry2d &self) { - bp::list leftpoints, rightpoints, leftdom, rightdom; + py::list leftpoints, rightpoints, leftdom, rightdom; for (int i = 0; i < self.splines.Size(); i++) { @@ -226,13 +227,13 @@ DLL_HEADER void ExportGeom2d() normal(0) = normal(1); normal(1) = -temp; - leftdom.append(self.GetSpline(i).leftdom); - rightdom.append(self.GetSpline(i).rightdom); + leftdom.append(py::cast(self.GetSpline(i).leftdom)); + rightdom.append(py::cast(self.GetSpline(i).rightdom)); - rightpoints.append(bp::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0)); - leftpoints.append(bp::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0)); + rightpoints.append(py::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0)); + leftpoints.append(py::make_tuple(point(0), point(1), normal(0)<0, normal(1)<0)); } - return bp::tuple(bp::make_tuple(leftpoints, rightpoints, leftdom, rightdom)); + return py::tuple(py::make_tuple(leftpoints, rightpoints, leftdom, rightdom)); })) .def("Print", FunctionPointer([](SplineGeometry2d &self) @@ -265,8 +266,10 @@ DLL_HEADER void ExportGeom2d() } -BOOST_PYTHON_MODULE(libgeom2d) { - ExportGeom2d(); +PYBIND11_PLUGIN(libgeom2d) { + py::module m("geom2d", "pybind geom2d"); + ExportGeom2d(m); + return m.ptr(); } #endif diff --git a/libsrc/gprim/geomfuncs.hpp b/libsrc/gprim/geomfuncs.hpp index 66cbca81..ea0070df 100644 --- a/libsrc/gprim/geomfuncs.hpp +++ b/libsrc/gprim/geomfuncs.hpp @@ -82,9 +82,10 @@ namespace netgen */ // inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2) - inline Vec<3> Cross (Vec<3> v1, Vec<3> v2) + template + inline Vec<3,T> Cross (Vec<3,T> v1, Vec<3,T> v2) { - return Vec<3> + return Vec<3,T> ( v1(1) * v2(2) - v1(2) * v2(1), v1(2) * v2(0) - v1(0) * v2(2), v1(0) * v2(1) - v1(1) * v2(0) ); diff --git a/libsrc/gprim/geomobjects.hpp b/libsrc/gprim/geomobjects.hpp index d4ca640f..e1e56785 100644 --- a/libsrc/gprim/geomobjects.hpp +++ b/libsrc/gprim/geomobjects.hpp @@ -17,7 +17,7 @@ namespace netgen template - class Point + class Point : public AlignedAlloc> { protected: @@ -39,16 +39,18 @@ namespace netgen Point (T ax, T ay, T az, T au) { x[0] = ax; x[1] = ay; x[2] = az; x[3] = au;} - Point (const Point & p2) - { for (int i = 0; i < D; i++) x[i] = p2.x[i]; } + template + Point (const Point & p2) + { for (int i = 0; i < D; i++) x[i] = p2(i); } explicit Point (const Vec & v) { for (int i = 0; i < D; i++) x[i] = v(i); } - Point & operator= (const Point & p2) + template + Point & operator= (const Point & p2) { - for (int i = 0; i < D; i++) x[i] = p2.x[i]; + for (int i = 0; i < D; i++) x[i] = p2(i); return *this; } @@ -65,7 +67,7 @@ namespace netgen }; template - class Vec + class Vec : public AlignedAlloc> { protected: @@ -90,17 +92,17 @@ namespace netgen Vec (const Vec & p2) { for (int i = 0; i < D; i++) x[i] = p2.x[i]; } - explicit Vec (const Point & p) + explicit Vec (const Point & p) { for (int i = 0; i < D; i++) x[i] = p(i); } - Vec (const Vec & p1, const Vec & p2) + Vec (const Vec & p1, const Vec & p2) { for(int i=0; i & p2) + template + Vec & operator= (const Vec & p2) { - for (int i = 0; i < D; i++) x[i] = p2.x[i]; + for (int i = 0; i < D; i++) x[i] = p2(i); return *this; } @@ -131,12 +133,12 @@ namespace netgen return l; } - const Vec & Normalize () + Vec & Normalize () { T l = Length(); - if (l != 0) - for (int i = 0; i < D; i++) - x[i] /= l; + // if (l != 0) + for (int i = 0; i < D; i++) + x[i] /= (l+1e-40); return *this; } @@ -148,7 +150,7 @@ namespace netgen template - class Mat + class Mat : public AlignedAlloc> { protected: diff --git a/libsrc/gprim/geomops.hpp b/libsrc/gprim/geomops.hpp index d08eef30..217c539b 100644 --- a/libsrc/gprim/geomops.hpp +++ b/libsrc/gprim/geomops.hpp @@ -29,10 +29,10 @@ namespace netgen - template - inline Point operator+ (const Point & a, const Vec & b) + template + inline Point operator+ (const Point & a, const Vec & b) { - Point res; + Point res; for (int i = 0; i < D; i++) res(i) = a(i) + b(i); return res; @@ -40,10 +40,10 @@ namespace netgen - template - inline Vec operator- (const Point & a, const Point & b) + template + inline Vec operator- (const Point & a, const Point & b) { - Vec res; + Vec res; for (int i = 0; i < D; i++) res(i) = a(i) - b(i); return res; @@ -69,10 +69,10 @@ namespace netgen - template - inline Vec operator* (double s, const Vec & b) + template + inline Vec operator* (T s, const Vec & b) { - Vec res; + Vec res; for (int i = 0; i < D; i++) res(i) = s * b(i); return res; diff --git a/libsrc/gprim/geomtest3d.hpp b/libsrc/gprim/geomtest3d.hpp index 05387da0..a8bb0266 100644 --- a/libsrc/gprim/geomtest3d.hpp +++ b/libsrc/gprim/geomtest3d.hpp @@ -78,6 +78,15 @@ extern double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3 extern double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, const Point3d & tp3, const Point3d & p); + inline double MinDistTP2 (const Point<2> & tp1, const Point<2> & tp2, + const Point<2> & tp3, const Point<2> & p) + { + return MinDistTP2 (Point<3> (tp1(0), tp1(1),0), + Point<3> (tp2(0), tp2(1),0), + Point<3> (tp3(0), tp3(1),0), + Point<3> (p(0), p(1),0)); + } + /// Minimal distance of the 2 lines [l1p1,l1p2] and [l2p1,l2p2] extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, const Point3d & l2p1, const Point3d & l2p2); diff --git a/libsrc/gprim/spline.cpp b/libsrc/gprim/spline.cpp index cc4c44d7..9736a499 100644 --- a/libsrc/gprim/spline.cpp +++ b/libsrc/gprim/spline.cpp @@ -189,6 +189,51 @@ namespace netgen if (tang * gradn < 0) u *= -1; } + + template + void SplineSeg3 :: GetCoeff (Vector & u, Point pref) const + { + DenseMatrix a(6, 6); + DenseMatrix ata(6, 6); + Vector f(6); + + u.SetSize(6); + + // ata.SetSymmetric(1); + + double t = 0; + for (int i = 0; i < 5; i++, t += 0.25) + { + Vec p = GetPoint (t)-pref; + a(i, 0) = p(0) * p(0); + a(i, 1) = p(1) * p(1); + a(i, 2) = p(0) * p(1); + a(i, 3) = p(0); + a(i, 4) = p(1); + a(i, 5) = 1; + } + a(5, 0) = 1; + + CalcAtA (a, ata); + + u = 0; + u(5) = 1; + a.MultTrans (u, f); + ata.Solve (f, u); + + // the sign + Point p0 = GetPoint(0); + Vec ht = GetTangent(0); + Vec<2> tang(ht(0), ht(1)); + + double gradx = u(3); + double grady = u(4); + // double gradx = 2.*u(0)*p0(0) + u(2)*p0(1) + u(3); + // double grady = 2.*u(1)*p0(1) + u(2)*p0(0) + u(4); + Vec<2> gradn (grady, -gradx); + + if (tang * gradn < 0) u *= -1; + } @@ -197,11 +242,25 @@ namespace netgen { double t_old = -1; + /* if(proj_latest_t > 0. && proj_latest_t < 1.) t = proj_latest_t; else t = 0.5; - + */ + double tmin = 1; + double dist_min2 = Dist2 (GetPoint(tmin), point); + for (double ti = 0; ti < 0.99; ti += 0.25) + { + double di = Dist2(GetPoint(ti), point); + if (di < dist_min2) + { + tmin = ti; + dist_min2 = di; + } + } + t = tmin; + Point phi; Vec phip,phipp,phimp; diff --git a/libsrc/gprim/spline.hpp b/libsrc/gprim/spline.hpp index a872d6e9..4121fd61 100644 --- a/libsrc/gprim/spline.hpp +++ b/libsrc/gprim/spline.hpp @@ -85,6 +85,7 @@ namespace netgen void PrintCoeff (ostream & ost) const; virtual void GetCoeff (Vector & coeffs) const = 0; + virtual void GetCoeff (Vector & coeffs, Point p0) const { ; } virtual void GetPoints (int n, Array > & points) const; @@ -95,6 +96,9 @@ namespace netgen Array < Point > & points, const double eps) const {points.SetSize(0);} + // is the point in the convex hull (increased by eps) of the spline ? + virtual bool InConvexHull (Point p, double eps) const = 0; + virtual double MaxCurvature(void) const = 0; virtual string GetType(void) const {return "splinebase";} @@ -135,12 +139,18 @@ namespace netgen virtual const GeomPoint & EndPI () const { return p2; } /// virtual void GetCoeff (Vector & coeffs) const; - + virtual void GetCoeff (Vector & coeffs, Point p0) const; + virtual string GetType(void) const {return "line";} virtual void LineIntersections (const double a, const double b, const double c, Array < Point > & points, const double eps) const; - + + virtual bool InConvexHull (Point p, double eps) const + { + return MinDistLP2 (p1, p2, p) < sqr(eps); + } + virtual double MaxCurvature(void) const {return 0;} virtual void Project (const Point point, Point & point_on_curve, double & t) const; @@ -178,7 +188,8 @@ namespace netgen virtual const GeomPoint & EndPI () const { return p3; } /// virtual void GetCoeff (Vector & coeffs) const; - + virtual void GetCoeff (Vector & coeffs, Point p0) const; + virtual string GetType(void) const {return "spline3";} const GeomPoint & TangentPoint (void) const { return p2; } @@ -186,6 +197,11 @@ namespace netgen DLL_HEADER virtual void LineIntersections (const double a, const double b, const double c, Array < Point > & points, const double eps) const; + virtual bool InConvexHull (Point p, double eps) const + { + return MinDistTP2 (p1, p2, p3, p) < sqr(eps); + } + DLL_HEADER virtual double MaxCurvature(void) const; DLL_HEADER virtual void Project (const Point point, Point & point_on_curve, double & t) const; @@ -232,6 +248,11 @@ namespace netgen virtual void LineIntersections (const double a, const double b, const double c, Array < Point > & points, const double eps) const; + virtual bool InConvexHull (Point p, double eps) const + { + return (Dist2 (p, pm) < sqr(radius+eps)); + } + virtual double MaxCurvature(void) const {return 1./radius;} }; @@ -261,6 +282,10 @@ namespace netgen virtual void GetCoeff (Vector & coeffs) const {;} virtual double MaxCurvature(void) const {return 1;} + + // needs implementation ... + virtual bool InConvexHull (Point p, double eps) const + { return true; } }; @@ -372,6 +397,19 @@ namespace netgen coeffs[5] = -dx * p1(1) + dy * p1(0); } + template + void LineSeg :: GetCoeff (Vector & coeffs, Point p) const + { + coeffs.SetSize(6); + + double dx = p2(0) - p1(0); + double dy = p2(1) - p1(1); + + coeffs[0] = coeffs[1] = coeffs[2] = 0; + coeffs[3] = -dy; + coeffs[4] = dx; + coeffs[5] = -dx * (p1(1)-p(1)) + dy * (p1(0)-p(0)); + } template @@ -597,6 +635,10 @@ namespace netgen virtual void GetCoeff (Vector & coeffs) const {;} virtual double MaxCurvature(void) const {return 1;} + + // needs implementation ... + virtual bool InConvexHull (Point p, double eps) const + { return true; } }; // Constructor diff --git a/libsrc/gprim/transform3d.hpp b/libsrc/gprim/transform3d.hpp index cd412433..49ce29ab 100644 --- a/libsrc/gprim/transform3d.hpp +++ b/libsrc/gprim/transform3d.hpp @@ -182,6 +182,9 @@ public: { to = m * from; } + + Point operator() (Point from) const { Point to; Transform(from, to); return to; } + Vec operator() (Vec from) const { Vec to; Transform(from, to); return to; } }; template diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index b7689b8a..7ffb837d 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -64,6 +64,10 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const Ng_Element ret; ret.type = NG_ELEMENT_TYPE(el.GetType()); ret.index = el.si; + if (mesh->GetDimension() == 2) + ret.mat = mesh->GetBCNamePtr(el.si-1); + else + ret.mat = nullptr; ret.points.num = el.GetNP(); ret.points.ptr = (int*)&(el[0]); diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index cd694c22..6f05ed86 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -4176,6 +4176,7 @@ namespace netgen T * dxdxi, size_t sdxdxi) { // multipointtrafovar++; + /* static int timer = NgProfiler::CreateTimer ("calcmultipointelementtrafo"); static int timer1 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 1"); static int timer2 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 2"); @@ -4183,6 +4184,7 @@ namespace netgen static int timer4 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 4"); static int timer5 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 5"); NgProfiler::RegionTimer reg(timer); + */ // NgProfiler::StartTimer (timer); // NgProfiler::StartTimer (timer1); if (mesh.coarsemesh) diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 3687e240..522eb8c8 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -1020,13 +1020,6 @@ namespace netgen class DLL_HEADER MeshingParameters { public: - - - - - - - /** 3d optimization strategy: // m .. move nodes @@ -1130,12 +1123,28 @@ namespace netgen MeshingParameters (); /// MeshingParameters (const MeshingParameters & mp2) = default; + MeshingParameters (MeshingParameters && mp2) = default; + MeshingParameters & operator= (const MeshingParameters & mp2) = default; + MeshingParameters & operator= (MeshingParameters && mp2) = default; /// void Print (ostream & ost) const; /// // void CopyFrom(const MeshingParameters & other); - + class MeshSizePoint + { + public: + Point<3> pnt; + double h; + MeshSizePoint (Point<3> _pnt, double _h) : pnt(_pnt), h(_h) { ; } + MeshSizePoint () = default; + MeshSizePoint (const MeshSizePoint &) = default; + MeshSizePoint (MeshSizePoint &&) = default; + MeshSizePoint & operator= (const MeshSizePoint &) = default; + MeshSizePoint & operator= (MeshSizePoint &&) = default; + }; + Array meshsize_points; + void (*render_function)(bool) = NULL; void Render(bool blocking = false) { diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index f8520800..cf99678f 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -1,7 +1,5 @@ #ifdef NG_PYTHON -#include -#include #include <../general/ngpython.hpp> #include @@ -9,10 +7,11 @@ #include #include #include <../interface/writeuser.hpp> +#include +#include using namespace netgen; -namespace bp = boost::python; namespace netgen { @@ -21,23 +20,19 @@ namespace netgen template -void ExportArray () +void ExportArray (py::module &m) { string name = string("Array_") + typeid(T).name(); - bp::class_,boost::noncopyable>(name.c_str()) - .def ("__len__", &Array::Size) + py::class_>(m, name.c_str()) + .def ("__len__", [] ( Array &self ) { return self.Size(); } ) .def ("__getitem__", FunctionPointer ([](Array & self, TIND i) -> T& { if (i < BASE || i >= BASE+self.Size()) - bp::exec("raise IndexError()\n"); + throw py::index_error(); return self[i]; }), - bp::return_value_policy()) - - .def ("__iter__", - bp::range (FunctionPointer([](Array & self) { return &self[BASE]; }), - FunctionPointer([](Array & self) { return &self[BASE+self.Size()]; }))) + py::return_value_policy::reference) ; } @@ -49,45 +44,46 @@ void TranslateException (const NgException & ex) } -DLL_HEADER void ExportNetgenMeshing() +DLL_HEADER void ExportNetgenMeshing(py::module &m) { - ModuleScope module("meshing"); - - bp::register_exception_translator(&TranslateException); - bp::class_("PointId", bp::init()) + py::class_(m, "PointId") + .def(py::init()) .def("__repr__", &ToString) .def("__str__", &ToString) - .add_property("nr", &PointIndex::operator int) + .def_property_readonly("nr", &PointIndex::operator int) .def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other) { return static_cast(self)==static_cast(other); }) ) .def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast(self); }) ) ; - bp::class_("ElementId3D", bp::init()) + py::class_(m, "ElementId3D") + .def(py::init()) .def("__repr__", &ToString) .def("__str__", &ToString) - .add_property("nr", &ElementIndex::operator int) + .def_property_readonly("nr", &ElementIndex::operator int) .def("__eq__" , FunctionPointer( [](ElementIndex &self, ElementIndex &other) { return static_cast(self)==static_cast(other); }) ) .def("__hash__" , FunctionPointer( [](ElementIndex &self ) { return static_cast(self); }) ) ; - bp::class_("ElementId2D", bp::init()) + py::class_(m, "ElementId2D") + .def(py::init()) .def("__repr__", &ToString) .def("__str__", &ToString) - .add_property("nr", &SurfaceElementIndex::operator int) + .def_property_readonly("nr", &SurfaceElementIndex::operator int) .def("__eq__" , FunctionPointer( [](SurfaceElementIndex &self, SurfaceElementIndex &other) { return static_cast(self)==static_cast(other); }) ) .def("__hash__" , FunctionPointer( [](SurfaceElementIndex &self ) { return static_cast(self); }) ) ; - bp::class_("ElementId1D", bp::init()) + py::class_(m, "ElementId1D") + .def(py::init()) .def("__repr__", &ToString) .def("__str__", &ToString) - .add_property("nr", &SegmentIndex::operator int) + .def_property_readonly("nr", &SegmentIndex::operator int) .def("__eq__" , FunctionPointer( [](SegmentIndex &self, SegmentIndex &other) { return static_cast(self)==static_cast(other); }) ) .def("__hash__" , FunctionPointer( [](SegmentIndex &self ) { return static_cast(self); }) ) @@ -96,176 +92,162 @@ DLL_HEADER void ExportNetgenMeshing() /* - bp::class_> ("Point") - .def(bp::init()) + py::class_> ("Point") + .def(py::init()) ; */ - bp::class_> */ >("MeshPoint") - .def(bp::init>()) + py::class_> */ >(m, "MeshPoint") + .def(py::init>()) .def("__str__", &ToString) .def("__repr__", &ToString) - .add_property("p", FunctionPointer([](const MeshPoint & self) - { - bp::list l; - l.append ( self[0] ); - l.append ( self[1] ); - l.append ( self[2] ); - return bp::tuple(l); + .def_property_readonly("p", FunctionPointer([](const MeshPoint & self) + { + py::list l; + l.append ( py::cast(self[0]) ); + l.append ( py::cast(self[1]) ); + l.append ( py::cast(self[2]) ); + return py::tuple(l); })) .def("__getitem__", FunctionPointer([](const MeshPoint & self, int index) { if(index<0 || index>2) - bp::exec("raise IndexError()\n"); + throw py::index_error(); return self[index]; })) ; - bp::class_("Element3D") - .def("__init__", bp::make_constructor - (FunctionPointer ([](int index, bp::list vertices) + py::class_(m, "Element3D") + .def("__init__", [](Element *instance, int index, py::list vertices) { - if (bp::len(vertices) == 4) + if (py::len(vertices) == 4) { - Element * tmp = new Element(TET); + new (instance) Element(TET); for (int i = 0; i < 4; i++) - (*tmp)[i] = bp::extract(vertices[i]); - tmp->SetIndex(index); - return tmp; + (*instance)[i] = py::extract(vertices[i])(); + instance->SetIndex(index); } - if (bp::len(vertices) == 6) + if (py::len(vertices) == 6) { - Element * tmp = new Element(PRISM); + new (instance) Element(PRISM); for (int i = 0; i < 6; i++) - (*tmp)[i] = bp::extract(vertices[i]); - tmp->SetIndex(index); - return tmp; + (*instance)[i] = py::extract(vertices[i])(); + instance->SetIndex(index); } - if (bp::len(vertices) == 8) + if (py::len(vertices) == 8) { - Element * tmp = new Element(HEX); + new (instance) Element(HEX); for (int i = 0; i < 8; i++) - (*tmp)[i] = bp::extract(vertices[i]); - tmp->SetIndex(index); - return tmp; + (*instance)[i] = py::extract(vertices[i])(); + instance->SetIndex(index); } throw NgException ("cannot create element"); - }), - bp::default_call_policies(), // need it to use arguments - (bp::arg("index")=1,bp::arg("vertices"))), + }, + py::arg("index")=1,py::arg("vertices"), "create volume element" ) .def("__repr__", &ToString) - .add_property("index", &Element::GetIndex, &Element::SetIndex) - .add_property("vertices", - FunctionPointer ([](const Element & self) -> bp::list + .def_property("index", &Element::GetIndex, &Element::SetIndex) + .def_property_readonly("vertices", + FunctionPointer ([](const Element & self) -> py::list { - bp::list li; + py::list li; for (int i = 0; i < self.GetNV(); i++) - li.append (self[i]); + li.append (py::cast(self[i])); return li; })) ; - bp::class_("Element2D") - .def("__init__", bp::make_constructor - (FunctionPointer ([](int index, bp::list vertices) + py::class_(m, "Element2D") + .def("__init__", + [](Element2d *instance, int index, py::list vertices) { - if (bp::len(vertices) == 3) + if (py::len(vertices) == 3) { - Element2d * tmp = new Element2d(TRIG); + new (instance) Element2d(TRIG); for (int i = 0; i < 3; i++) - (*tmp)[i] = bp::extract(vertices[i]); - tmp->SetIndex(index); - return tmp; + (*instance)[i] = py::extract(vertices[i])(); + instance->SetIndex(index); } else { - Element2d * tmp = new Element2d(QUAD); + new (instance) Element2d(QUAD); for (int i = 0; i < 4; i++) - (*tmp)[i] = bp::extract(vertices[i]); - tmp->SetIndex(index); - return tmp; + (*instance)[i] = py::extract(vertices[i])(); + instance->SetIndex(index); } - }), - bp::default_call_policies(), // need it to use arguments - (bp::arg("index")=1,bp::arg("vertices"))), + }, + py::arg("index")=1,py::arg("vertices"), "create surface element" ) - .add_property("index", &Element2d::GetIndex, &Element2d::SetIndex) - .add_property("vertices", - FunctionPointer([](const Element2d & self) -> bp::list + .def_property("index", &Element2d::GetIndex, &Element2d::SetIndex) + .def_property_readonly("vertices", + FunctionPointer([](const Element2d & self) -> py::list { - bp::list li; + py::list li; for (int i = 0; i < self.GetNV(); i++) - li.append(self[i]); + li.append(py::cast(self[i])); return li; })) ; - bp::class_("Element1D") - .def("__init__", bp::make_constructor - (FunctionPointer ([](bp::list vertices, bp::list surfaces, int index) + py::class_(m, "Element1D") + .def("__init__", + [](Segment *instance, py::list vertices, py::list surfaces, int index) { - Segment * tmp = new Segment; + new (instance) Segment(); for (int i = 0; i < 2; i++) - (*tmp)[i] = bp::extract(vertices[i]); - tmp -> si = index; + (*instance)[i] = py::extract(vertices[i])(); + instance -> si = index; if (len(surfaces)) { - tmp->surfnr1 = bp::extract(surfaces[0]); - tmp->surfnr2 = bp::extract(surfaces[1]); + instance->surfnr1 = py::extract(surfaces[0])(); + instance->surfnr2 = py::extract(surfaces[1])(); } - return tmp; - }), - bp::default_call_policies(), - (bp::arg("vertices"), - bp::arg("surfaces")=bp::list(), - bp::arg("index")=1 - )), + }, + py::arg("vertices"), + py::arg("surfaces")=py::list(), + py::arg("index")=1, "create segment element" ) .def("__repr__", &ToString) - .add_property("vertices", - FunctionPointer ([](const Segment & self) -> bp::list + .def_property_readonly("vertices", + FunctionPointer ([](const Segment & self) -> py::list { - bp::list li; + py::list li; for (int i = 0; i < 2; i++) - li.append (self[i]); + li.append (py::cast(self[i])); return li; })) - .add_property("surfaces", - FunctionPointer ([](const Segment & self) -> bp::list + .def_property_readonly("surfaces", + FunctionPointer ([](const Segment & self) -> py::list { - bp::list li; - li.append (self.surfnr1); - li.append (self.surfnr2); + py::list li; + li.append (py::cast(self.surfnr1)); + li.append (py::cast(self.surfnr2)); return li; })) ; - bp::class_("Element0D") - .def("__init__", bp::make_constructor - (FunctionPointer ([](PointIndex vertex, int index) + py::class_(m, "Element0D") + .def("__init__", + [](Element0d *instance, PointIndex vertex, int index) { - Element0d * tmp = new Element0d; - tmp->pnum = vertex; - tmp->index = index; - return tmp; - }), - bp::default_call_policies(), - (bp::arg("vertex"), - bp::arg("index")=1 - )), + new (instance) Element0d; + instance->pnum = vertex; + instance->index = index; + }, + py::arg("vertex"), + py::arg("index")=1, "create point element" ) .def("__repr__", &ToString) - .add_property("vertices", - FunctionPointer ([](const Element0d & self) -> bp::list + .def_property_readonly("vertices", + FunctionPointer ([](const Element0d & self) -> py::list { - bp::list li; - li.append (self.pnum); + py::list li; + li.append (py::cast(self.pnum)); return li; })) ; @@ -274,62 +256,65 @@ DLL_HEADER void ExportNetgenMeshing() - bp::class_("FaceDescriptor") - .def(bp::init()) - .def("__init__", bp::make_constructor - (FunctionPointer ([](int surfnr, int domin, int domout, int bc) + py::class_(m, "FaceDescriptor") + .def(py::init()) + .def("__init__", + [](FaceDescriptor *instance, int surfnr, int domin, int domout, int bc) { - auto fd = new FaceDescriptor(); - fd->SetSurfNr(surfnr); - fd->SetDomainIn(domin); - fd->SetDomainOut(domout); - fd->SetBCProperty(bc); - return fd; - }), - bp::default_call_policies(), // need it to use arguments - (bp::arg("surfnr")=1, - bp::arg("domin")=1, - bp::arg("domout")=0, - bp::arg("bc")=0 - )), + new (instance) FaceDescriptor(); + instance->SetSurfNr(surfnr); + instance->SetDomainIn(domin); + instance->SetDomainOut(domout); + instance->SetBCProperty(bc); + }, + py::arg("surfnr")=1, + py::arg("domin")=1, + py::arg("domout")=py::int_(0), + py::arg("bc")=py::int_(0), "create facedescriptor") .def("__str__", &ToString) .def("__repr__", &ToString) - .add_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr) - .add_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn) - .add_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut) - .add_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty) - .add_property("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); })) + .def_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr) + .def_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn) + .def_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut) + .def_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty) + .def_property_readonly("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); })) + .def("SetSurfaceColor", [](FaceDescriptor & self, py::list color ) + { + Vec3d c; + c.X() = py::extract(color[0])(); + c.Y() = py::extract(color[1])(); + c.Z() = py::extract(color[2])(); + self.SetSurfColour(c); + }) ; - ExportArray(); - ExportArray(); - ExportArray(); - ExportArray(); - ExportArray(); - ExportArray(); - ; - -#if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100) - bp::register_ptr_to_python>(); -#endif - bp::class_,boost::noncopyable>("Mesh", bp::no_init) - // .def(bp::init<>("create empty mesh")) + ExportArray(m); + ExportArray(m); + ExportArray(m); + ExportArray(m); + ExportArray(m); + ExportArray(m); - .def("__init__", bp::make_constructor - (FunctionPointer ([](int dim) + py::implicitly_convertible< int, PointIndex>(); + + py::class_(m, "Ngx_Mesh") + .def(py::init>()) + .def_property_readonly("ngmesh", &Ngx_Mesh::GetMesh) + ; + py::class_>(m, "Mesh") + // .def(py::init<>("create empty mesh")) + + .def("__init__", + [](Mesh *instance, int dim) { - auto mesh = make_shared(); - mesh->SetDimension(dim); - return mesh; - }), - bp::default_call_policies(), // need it to use named arguments - ( - bp::arg("dim")=3 - ) - )) + new (instance) Mesh(); + instance->SetDimension(dim); + }, + py::arg("dim")=3 + ) .def("__str__", &ToString) @@ -371,8 +356,8 @@ DLL_HEADER void ExportNetgenMeshing() })) // static_cast(&Mesh::Load)) .def("Save", static_cast(&Mesh::Save)) - .def("Export", FunctionPointer - ([] (Mesh & self, string filename, string format) + .def("Export", + [] (Mesh & self, string filename, string format) { if (WriteUserFormat (format, self, *self.GetGeometry(), filename)) { @@ -384,35 +369,35 @@ DLL_HEADER void ExportNetgenMeshing() err += string("'") + name + "'\n"; throw NgException (err); } - }), - (bp::arg("self"), bp::arg("filename"), bp::arg("format"))) + }, + py::arg("filename"), py::arg("format")) - .add_property("dim", &Mesh::GetDimension, &Mesh::SetDimension) + .def_property("dim", &Mesh::GetDimension, &Mesh::SetDimension) .def("Elements3D", static_cast&(Mesh::*)()> (&Mesh::VolumeElements), - bp::return_value_policy()) + py::return_value_policy::reference) .def("Elements2D", static_cast&(Mesh::*)()> (&Mesh::SurfaceElements), - bp::return_value_policy()) + py::return_value_policy::reference) .def("Elements1D", static_cast&(Mesh::*)()> (&Mesh::LineSegments), - bp::return_value_policy()) + py::return_value_policy::reference) .def("Elements0D", FunctionPointer([] (Mesh & self) -> Array& { return self.pointelements; } ), - bp::return_value_policy()) + py::return_value_policy::reference) .def("Points", static_cast (&Mesh::Points), - bp::return_value_policy()) + py::return_value_policy::reference) .def("FaceDescriptor", static_cast (&Mesh::GetFaceDescriptor), - bp::return_value_policy()) + py::return_value_policy::reference) .def("GetNFaceDescriptors", &Mesh::GetNFD) @@ -457,16 +442,23 @@ DLL_HEADER void ExportNetgenMeshing() .def ("SetMaterial", &Mesh::SetMaterial) .def ("GetMaterial", FunctionPointer([](Mesh & self, int domnr) { return string(self.GetMaterial(domnr)); })) - - .def ("GenerateVolumeMesh", FunctionPointer - ([](Mesh & self) + + .def ("GenerateVolumeMesh", + [](Mesh & self, py::object pymp) { cout << "generate vol mesh" << endl; + MeshingParameters mp; - mp.optsteps3d = 5; + if (py::extract(pymp).check()) + mp = py::extract(pymp)(); + else + { + mp.optsteps3d = 5; + } MeshVolume (mp, self); OptimizeVolume (mp, self); - })) + }, + py::arg("mp")=NGDummyArgument()) .def ("OptimizeVolumeMesh", FunctionPointer ([](Mesh & self) @@ -491,18 +483,19 @@ DLL_HEADER void ExportNetgenMeshing() self.SetGeometry(geo); })) - .def ("SetGeometry", FunctionPointer - ([](Mesh & self, shared_ptr geo) - { - self.SetGeometry(geo); - })) + // TODO: fix this dependency on libgeom2d.so +// .def ("SetGeometry", FunctionPointer +// ([](Mesh & self, shared_ptr geo) +// { +// self.SetGeometry(geo); +// })) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree) .def ("BoundaryLayer", FunctionPointer - ([](Mesh & self, int bc, bp::list thicknesses, int volnr, bp::list materials) + ([](Mesh & self, int bc, py::list thicknesses, int volnr, py::list materials) { - int n = bp::len(thicknesses); + int n = py::len(thicknesses); BoundaryLayerParameters blp; for (int i = 1; i <= self.GetNFD(); i++) @@ -521,9 +514,9 @@ DLL_HEADER void ExportNetgenMeshing() cout << "maxind = " << maxind << endl; for ( int i=0; i(thicknesses[i])()) ; + blp.heights.Append( py::extract(thicknesses[i])()) ; blp.new_matnrs.Append( maxind+1+i ); - self.SetMaterial (maxind+1+i, bp::extract(materials[i])().c_str()); + self.SetMaterial (maxind+1+i, py::extract(materials[i])().c_str()); } blp.bulk_matnr = volnr; GenerateBoundaryLayer (self, blp); @@ -558,39 +551,48 @@ DLL_HEADER void ExportNetgenMeshing() )) ; + py::implicitly_convertible< shared_ptr, Ngx_Mesh >(); typedef MeshingParameters MP; - bp::class_ ("MeshingParameters", bp::init<>()) - .def("__init__", bp::make_constructor - (FunctionPointer ([](double maxh, bool quad_dominated) + py::class_ (m, "MeshingParameters") + .def(py::init<>()) + .def("__init__", + [](MP *instance, double maxh, bool quad_dominated, int optsteps2d, int optsteps3d) { - auto tmp = new MeshingParameters; - tmp->maxh = maxh; - tmp->quad = int(quad_dominated); - return tmp; - }), - bp::default_call_policies(), // need it to use arguments - ( - bp::arg("maxh")=1000, - bp::arg("quad_dominated")=false - )), + new (instance) MeshingParameters; + instance->maxh = maxh; + instance->quad = int(quad_dominated); + instance->optsteps2d = optsteps2d; + instance->optsteps3d = optsteps3d; + }, + py::arg("maxh")=1000, + py::arg("quad_dominated")=false, + py::arg("optsteps2d") = 3, + py::arg("optsteps3d") = 3 + , "create meshing parameters" ) .def("__str__", &ToString) - .add_property("maxh", + .def_property("maxh", FunctionPointer ([](const MP & mp ) { return mp.maxh; }), FunctionPointer ([](MP & mp, double maxh) { return mp.maxh = maxh; })) - + .def("RestrictH", FunctionPointer + ([](MP & mp, double x, double y, double z, double h) + { + mp.meshsize_points.Append ( MeshingParameters::MeshSizePoint (Point<3> (x,y,z), h)); + }), + py::arg("x"), py::arg("y"), py::arg("z"), py::arg("h") + ) ; - bp::def("SetTestoutFile", FunctionPointer ([] (const string & filename) + m.def("SetTestoutFile", FunctionPointer ([] (const string & filename) { delete testout; testout = new ofstream (filename); })); - bp::def("SetMessageImportance", FunctionPointer ([] (int importance) + m.def("SetMessageImportance", FunctionPointer ([] (int importance) { int old = printmessage_importance; printmessage_importance = importance; @@ -598,14 +600,11 @@ DLL_HEADER void ExportNetgenMeshing() })); } - - -BOOST_PYTHON_MODULE(libmesh) { - ExportNetgenMeshing(); +PYBIND11_PLUGIN(libmesh) { + py::module m("mesh", "pybind mesh"); + ExportNetgenMeshing(m); + return m.ptr(); } - - - #endif diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index 8e122fd3..01346f40 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -62,7 +62,63 @@ namespace netgen delete vert2pointelement; } + template + void LoopOverEdges (const Mesh & mesh, MeshTopology & top, PointIndex v, + FUNC func) + { + for (ElementIndex elnr : top.GetVertexElements(v)) + { + const Element & el = mesh[elnr]; + int neledges = MeshTopology::GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + // edge.Sort(); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + if (edge.I1() != v) continue; + + func (edge, elnr, k, 3, edgedir); + } + } + + for (SurfaceElementIndex elnr : top.GetVertexSurfaceElements(v)) + { + const Element2d & el = mesh[elnr]; + + int neledges = MeshTopology::GetNEdges (el.GetType()); + const ELEMENT_EDGE * eledges = MeshTopology::GetEdges0 (el.GetType()); + + for (int k = 0; k < neledges; k++) + { + INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); + // edge.Sort(); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + if (edge.I1() != v) continue; + + func (edge, elnr, k, 2, edgedir); + } + } + + for (SegmentIndex elnr : top.GetVertexSegments(v)) + { + const Segment & el = mesh[elnr]; + INDEX_2 edge(el[0], el[1]); + int edgedir = (edge.I1() > edge.I2()); + if (edgedir) swap (edge.I1(), edge.I2()); + + edge.Sort(); + if (edge.I1() != v) continue; + + func (edge, elnr, 0, 1, edgedir); + } + } + template void LoopOverFaces (const Mesh & mesh, MeshTopology & top, PointIndex v, FUNC func) @@ -139,11 +195,12 @@ namespace netgen swap (face4.I2(), face4.I4()); } + if (face4.I1() != v) continue; + func(face4, elnr, j, true, facedir); /* INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); - if (face.I1() != v) continue; if (vert2face.Used (face)) { @@ -256,10 +313,10 @@ namespace netgen swap (face4.I2(), face4.I4()); } + if (face4.I1() != v) continue; func(face4, elnr, 0, false, facedir); /* INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); - if (face.I1() != v) continue; if (vert2face.Used (face)) facenum = vert2face.Get(face); @@ -300,7 +357,6 @@ namespace netgen int nseg = mesh.GetNSeg(); int np = mesh.GetNP(); int nv = mesh.GetNV(); - int ned = edge2vert.Size(); if (id == 0) PrintMessage (3, "Update mesh topology"); @@ -417,7 +473,7 @@ namespace netgen cnt[edge2vert[i][0]]++; TABLE vert2edge (cnt); for (int i = 0; i < edge2vert.Size(); i++) - vert2edge.AddSave (edge2vert[i][0], i+1); + vert2edge.AddSave (edge2vert[i][0], i); // ensure all coarse grid and intermediate level edges cnt = 0; @@ -433,167 +489,141 @@ namespace netgen if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]); } - Array edgenr(nv); - Array edgeflag(nv); - Array vertex2; - edgeflag = PointIndex::BASE-1; - - ned = edge2vert.Size(); + int max_edge_on_vertex = 0; for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) { - vertex2.SetSize (0); - - for (int j = 0; j < vert2edge[i].Size(); j++) - { - int ednr = vert2edge[i][j]; - int i2 = edge2vert.Get(ednr)[1]; - edgeflag[i2] = i; - edgenr[i2] = ednr; - } - - for (int j = 0; j < vert2vertcoarse[i].Size(); j++) - { - int v2 = vert2vertcoarse[i][j]; - if (edgeflag[v2] < i) - { - edgeflag[v2] = i; - vertex2.Append (v2); - } - } - - FlatArray v2els = (*vert2element)[i]; - for (int j = 0; j < v2els.Size(); j++) - { - const Element & el = mesh[v2els[j]]; - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - edge.Sort(); - if (edge.I1() != i) continue; - - if (edgeflag[edge.I2()] < i) - { - vertex2.Append (edge.I2()); - edgeflag[edge.I2()] = i; - } - } - } - - for (int j = 0; j < (*vert2surfelement)[i].Size(); j++) - { - SurfaceElementIndex elnr = (*vert2surfelement)[i][j]; - const Element2d & el = mesh.SurfaceElement (elnr); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - edge.Sort(); - if (edge.I1() != i) continue; - - if (edgeflag[edge.I2()] < i) - { - vertex2.Append (edge.I2()); - edgeflag[edge.I2()] = i; - } - } - } - - for (int j = 0; j < (*vert2segment)[i].Size(); j++) - { - SegmentIndex elnr = (*vert2segment)[i][j]; - const Segment & el = mesh.LineSegment (elnr); - - INDEX_2 edge(el[0], el[1]); - edge.Sort(); - if (edge.I1() != i) continue; - - if (edgeflag[edge.I2()] < i) - { - vertex2.Append (edge.I2()); - edgeflag[edge.I2()] = i; - } - } - - QuickSort (vertex2); - for (int j = 0; j < vertex2.Size(); j++) - { - edgenr[vertex2[j]] = ++ned; - edge2vert.Append (INDEX_2 (i, vertex2[j])); - } - - - for (int j = 0; j < (*vert2element)[i].Size(); j++) - { - ElementIndex elnr = (*vert2element)[i][j]; - const Element & el = mesh[elnr]; - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) continue; - - int edgenum = edgenr[edge.I2()]; - edges[elnr][k].nr = edgenum-1; - edges[elnr][k].orient = edgedir; - } - } - - for (int j = 0; j < (*vert2surfelement)[i].Size(); j++) - { - SurfaceElementIndex elnr = (*vert2surfelement)[i][j]; - const Element2d & el = mesh.SurfaceElement (elnr); - - int neledges = GetNEdges (el.GetType()); - const ELEMENT_EDGE * eledges = GetEdges0 (el.GetType()); - - for (int k = 0; k < neledges; k++) - { - INDEX_2 edge(el[eledges[k][0]], el[eledges[k][1]]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) continue; - - int edgenum = edgenr[edge.I2()]; - surfedges[elnr][k].nr = edgenum-1; - surfedges[elnr][k].orient = edgedir; - } - } - - for (int j = 0; j < (*vert2segment)[i].Size(); j++) - { - SegmentIndex elnr = (*vert2segment)[i][j]; - const Segment & el = mesh.LineSegment (elnr); - - INDEX_2 edge(el[0], el[1]); - - int edgedir = (edge.I1() > edge.I2()); - if (edgedir) swap (edge.I1(), edge.I2()); - - if (edge.I1() != i) continue; - - int edgenum = edgenr[edge.I2()]; - segedges[elnr].nr = edgenum-1; - segedges[elnr].orient = edgedir; - } + int onv = vert2edge[i].Size() + vert2vertcoarse[i].Size() + + 4*(*vert2element)[i].Size() + 2*(*vert2surfelement)[i].Size() + (*vert2segment)[i].Size(); + max_edge_on_vertex = max (onv, max_edge_on_vertex); } - } + + // count edges associated with vertices + cnt = 0; + + ParallelForRange + (tm, mesh.Points().Size(), + [&] (size_t begin, size_t end) + { + INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); + for (PointIndex v = begin+PointIndex::BASE; + v < end+PointIndex::BASE; v++) + { + v2eht.DeleteData(); + for (int ednr : vert2edge[v]) + { + int v2 = edge2vert[ednr][1]; + v2eht.Set (v2, ednr); + } + + int cnti = 0; + + for (int v2 : vert2vertcoarse[v]) + if (!v2eht.Used(v2)) + { + cnti++; + v2eht.Set (v2, 33); // some value + } + + LoopOverEdges (mesh, *this, v, + [&] (INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) + { + if (!v2eht.Used (edge.I2())) + { + cnti++; + v2eht.Set (edge.I2(), 33); // something + } + }); + cnt[v] = cnti; + } + } ); + + // accumulate number of edges + int ned = edge2vert.Size(); + for (auto v : mesh.Points().Range()) + { + auto hv = cnt[v]; + cnt[v] = ned; + ned += hv; + } + edge2vert.SetSize(ned); + + + // INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); + // Array vertex2; + // for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) + + ParallelForRange + (tm, mesh.Points().Size(), + [&] (size_t begin, size_t end) + { + INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); + Array vertex2; + for (PointIndex v = begin+PointIndex::BASE; + v < end+PointIndex::BASE; v++) + { + int ned = cnt[v]; + v2eht.DeleteData(); + vertex2.SetSize (0); + + for (int ednr : vert2edge[v]) + { + int v2 = edge2vert[ednr][1]; + v2eht.Set (v2, ednr); + } + + for (int v2 : vert2vertcoarse[v]) + if (!v2eht.Used(v2)) + { + v2eht.Set (v2, 33); // some value + vertex2.Append (v2); + } + + LoopOverEdges (mesh, *this, v, + [&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) + { + if (!v2eht.Used(edge.I2())) + { + vertex2.Append (edge.I2()); + v2eht.Set (edge.I2(), 33); + } + }); + + QuickSort (vertex2); + + for (int j = 0; j < vertex2.Size(); j++) + { + v2eht.Set (vertex2[j], ned); + edge2vert[ned] = INDEX_2 (v, vertex2[j]); + ned++; + } + + LoopOverEdges (mesh, *this, v, + [&](INDEX_2 edge, int elnr, int loc_edge, int element_dim, int edgedir) + { + int edgenum = v2eht.Get(edge.I2()); + switch (element_dim) + { + case 3: + edges[elnr][loc_edge].nr = edgenum; + edges[elnr][loc_edge].orient = edgedir; + break; + case 2: + surfedges[elnr][loc_edge].nr = edgenum; + surfedges[elnr][loc_edge].orient = edgedir; + break; + case 1: + segedges[elnr].nr = edgenum; + segedges[elnr].orient = edgedir; + break; + } + }); + } + } ); + } + // generate faces @@ -638,7 +668,6 @@ namespace netgen - NgProfiler::StopTimer (timer2a); NgProfiler::StartTimer (timer2b); @@ -764,8 +793,9 @@ namespace netgen }); } }); - /* - int oldnfa = face2vert.Size(); + + /* + int oldnfa = face2vert.Size(); int nfa = oldnfa; INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); @@ -1010,9 +1040,8 @@ namespace netgen } } } -*/ - face2vert.SetAllocSize (nfa); + */ // *testout << "face2vert = " << endl << face2vert << endl; diff --git a/libsrc/stlgeom/stltool.cpp b/libsrc/stlgeom/stltool.cpp index 9028cd72..855cec42 100644 --- a/libsrc/stlgeom/stltool.cpp +++ b/libsrc/stlgeom/stltool.cpp @@ -615,7 +615,7 @@ STLChart :: STLChart(STLGeometry * ageometry) geometry = ageometry; - if ( (stlparam.usesearchtree == 1)) + if ( stlparam.usesearchtree == 1) searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1), geometry->GetBoundingBox().PMax() + Vec3d(1,1,1)); else @@ -1044,7 +1044,8 @@ int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> } - if (possible && 0) + // if (possible && 0) + if (false) for (i = 0; i <= divisions; i++) { diff --git a/libsrc/visualization/mvdraw.cpp b/libsrc/visualization/mvdraw.cpp index 06e8b7e8..04d98e6c 100644 --- a/libsrc/visualization/mvdraw.cpp +++ b/libsrc/visualization/mvdraw.cpp @@ -776,13 +776,16 @@ namespace netgen { string displname; - // Display * dpy = glXGetCurrentDisplay(); + Display * dpy = glXGetCurrentDisplay(); GLXDrawable drawable = glXGetCurrentDrawable(); GLXContext ctx = glXGetCurrentContext(); GLXContextID xid = glXGetContextIDEXT (ctx); displname = XDisplayName (0); + if( glXIsDirect ( dpy, ctx ) ) + cout << "WARNING: direct rendering enabled; this might break mpi-parallel netgen (especially if X-forwarding is used! (to disable, change -indirect to true in ng/drawing.tcl)" << endl; + /* cout << "Init Parallel GL" << endl; cout << "DisplayName = " << displname << endl; diff --git a/libsrc/visualization/soldata.hpp b/libsrc/visualization/soldata.hpp index 1a7dfd3e..20028762 100644 --- a/libsrc/visualization/soldata.hpp +++ b/libsrc/visualization/soldata.hpp @@ -92,6 +92,18 @@ namespace netgen return res; } +#ifdef __AVX__ + virtual bool GetMultiSurfValue (size_t selnr, size_t facetnr, size_t npts, + const __m256d * xref, + const __m256d * x, + const __m256d * dxdxref, + __m256d * values) + { + cerr << "GetMultiSurfVaue not overloaded" << endl; + return false; + } +#endif + virtual bool GetSegmentValue (int segnr, double xref, double * values) { return false; } diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index 989fb9a9..2a58e6e6 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -3529,14 +3529,11 @@ namespace netgen #ifdef NG_PYTHON -#include #include <../general/ngpython.hpp> -namespace bp = boost::python; -DLL_HEADER void ExportMeshVis() +DLL_HEADER void ExportMeshVis(py::module &m) { using namespace netgen; - ModuleScope myscope("meshvis"); vispar.drawcolorbar = true; vispar.drawnetgenlogo = true; vispar.drawcoordinatecross = true; @@ -3546,12 +3543,12 @@ DLL_HEADER void ExportMeshVis() vispar.drawtets = true; vispar.drawprisms = true; vispar.drawoutline = true; - bp::class_> - ("VisualSceneMesh", bp::no_init) + py::class_> + (m, "VisualSceneMesh") .def("Draw", &VisualSceneMesh::DrawScene) ; - bp::def("VS", FunctionPointer + m.def("VS", FunctionPointer ([](shared_ptr mesh) { auto vs = make_shared(); @@ -3560,15 +3557,23 @@ DLL_HEADER void ExportMeshVis() return vs; })); - bp::def("MouseMove", FunctionPointer + m.def("MouseMove", FunctionPointer ([](VisualSceneMesh &vsmesh, int oldx, int oldy, int newx, int newy, char mode) { vsmesh.MouseMove(oldx, oldy, newx, newy, mode); })); + m.def("SelectFace", FunctionPointer + ([] (int facenr) { + vsmesh.SetSelectedFace(facenr); + })); + m.def("GetGlobalMesh", FunctionPointer + ([] () { + return vsmesh.GetMesh(); + })); } -BOOST_PYTHON_MODULE(libvisual) -{ - ExportMeshVis(); -} +// BOOST_PYTHON_MODULE(libvisual) +// { +// ExportMeshVis(); +// } #endif diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index fbe8d732..6a02e54e 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -8,10 +8,7 @@ // #include #include - #include - - namespace netgen { @@ -441,6 +438,20 @@ namespace netgen // glEnable(GL_BLEND); glDisable(GL_BLEND); glCallList (surfellist); + +#ifdef USE_BUFFERS + static int timer = NgProfiler::CreateTimer ("Solution::drawing - DrawSurfaceElements VBO"); + NgProfiler::StartTimer(timer); + glEnableClientState(GL_VERTEX_ARRAY); + glEnableClientState(GL_NORMAL_ARRAY); + glEnableClientState(GL_TEXTURE_COORD_ARRAY); + glDrawElements(GL_TRIANGLES, surfel_vbo_size, GL_UNSIGNED_INT, 0); + glDisableClientState(GL_VERTEX_ARRAY); + glDisableClientState(GL_NORMAL_ARRAY); + glDisableClientState(GL_TEXTURE_COORD_ARRAY); + NgProfiler::StopTimer(timer); +#endif + glDisable(GL_BLEND); /* // transparent test ... @@ -1182,6 +1193,17 @@ namespace netgen shared_ptr mesh = GetMesh(); static int timer = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements"); + static int timerstart = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements start"); + static int timerloops = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements loops"); + static int timerlist = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements list"); + static int timerbuffer = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements buffer"); + static int timer1 = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1"); + static int timer1a = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1a"); + static int timer1b = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1b"); + static int timer1c = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 1c"); + static int timer2 = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 2"); + static int timer2a = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 2a"); + static int timer2b = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements 2b"); NgProfiler::RegionTimer reg (timer); @@ -1213,7 +1235,7 @@ namespace netgen } #endif - + NgProfiler::StartTimer(timerstart); if (surfellist) glDeleteLists (surfellist, 1); @@ -1256,10 +1278,58 @@ namespace netgen Array values(npt); Array mvalues(npt); + int sol_comp = (sol && sol->draw_surface) ? sol->components : 0; +#ifdef __AVX__ + Array> > simd_pref ( (npt+SIMD::Size()-1)/SIMD::Size() ); + Array> > simd_points ( (npt+SIMD::Size()-1)/SIMD::Size() ); + Array> > simd_dxdxis ( (npt+SIMD::Size()-1)/SIMD::Size() ); + Array> > simd_nvs( (npt+SIMD::Size()-1)/SIMD::Size() ); + Array> simd_values( (npt+SIMD::Size()-1)/SIMD::Size() * sol_comp); + + +#endif + + // Array> glob_pnts; + // Array> glob_nvs; + // Array glob_values; + if (sol && sol->draw_surface) mvalues.SetSize (npt * sol->components); - + Array > valuesc(npt); + +#ifdef USE_BUFFERS + if (has_surfel_vbo) + glDeleteBuffers (4, &surfel_vbo[0]); + glGenBuffers (4, &surfel_vbo[0]); + has_surfel_vbo = true; + glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[0]); + glBufferData (GL_ARRAY_BUFFER, + nse*npt*sizeof(Point<3,double>), + NULL, GL_STATIC_DRAW); + glVertexPointer(3, GL_DOUBLE, 0, 0); + // glEnableClientState(GL_VERTEX_ARRAY); + + glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[1]); + glBufferData (GL_ARRAY_BUFFER, + nse*npt*sizeof(Vec<3,double>), + NULL, GL_STATIC_DRAW); + // glEnableClientState(GL_NORMAL_ARRAY); + glNormalPointer(GL_DOUBLE, 0, 0); + + // glEnableClientState(GL_TEXTURE_COORD_ARRAY); + glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[2]); + glBufferData (GL_ARRAY_BUFFER, nse*npt*sizeof(double), NULL, GL_STATIC_DRAW); + glTexCoordPointer(1, GL_DOUBLE, 0, 0); + + glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]); + glBufferData(GL_ELEMENT_ARRAY_BUFFER, nse*npt*6*sizeof(int), NULL, GL_STATIC_DRAW); + surfel_vbo_size = 0; +#endif + + + NgProfiler::StopTimer(timerstart); + for (SurfaceElementIndex sei = 0; sei < nse; sei++) { const Element2d & el = (*mesh)[sei]; @@ -1413,7 +1483,38 @@ namespace netgen n = 1 << subdivisions; double invn = 1.0 / n; npt = (n+1)*(n+2)/2; + NgProfiler::StartTimer(timerloops); + size_t base_pi = 0; +#ifdef __AVX__ + for (int iy = 0, ii = 0; iy <= n; iy++) + for (int ix = 0; ix <= n-iy; ix++, ii++) + pref[ii] = Point<2> (ix*invn, iy*invn); + + constexpr size_t simd_size = SIMD::Size(); + size_t simd_npt = (npt+simd_size-1)/simd_size; + + for (size_t i = 0; i < simd_npt; i++) + { + simd_pref[i](0).SIMD_function ([&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](0) : 0; }, std::true_type()); + simd_pref[i](1).SIMD_function ([&] (size_t j) { size_t ii = i*simd_size+j; return (ii < npt) ? pref[ii](1) : 0; }, std::true_type()); + } +#endif + + Array ind_reftrig; + for (int iy = 0, ii = 0; iy < n; iy++,ii++) + for (int ix = 0; ix < n-iy; ix++, ii++) + { + int nv = (ix+iy+1 < n) ? 6 : 3; + int ind[] = { ii, ii+1, ii+n-iy+1, + ii+n-iy+1, ii+1, ii+n-iy+2 }; + for (int j = 0; j < nv; j++) + ind_reftrig.Append (ind[j]); + } + Array glob_ind; + glob_ind.SetSize(ind_reftrig.Size()); + + for(SurfaceElementIndex sei = 0; sei < nse; sei++) { const Element2d & el = (*mesh)[sei]; @@ -1436,6 +1537,90 @@ namespace netgen if ( el.GetType() == TRIG || el.GetType() == TRIG6 ) { + NgProfiler::StartTimer(timer1); +#ifdef __AVX_try_it_out__ + // NgProfiler::StartTimer(timer1a); + bool curved = curv.IsSurfaceElementCurved(sei); + + if (curved) + { + mesh->GetCurvedElements(). + CalcMultiPointSurfaceTransformation<3> (sei, simd_npt, + &simd_pref[0](0), 2, + &simd_points[0](0), 3, + &simd_dxdxis[0](0,0), 6); + + for (size_t ii = 0; ii < simd_npt; ii++) + simd_nvs[ii] = Cross (simd_dxdxis[ii].Col(0), simd_dxdxis[ii].Col(1)).Normalize(); + } + else + { + Point<3,SIMD> p1 = mesh->Point (el[0]); + Point<3,SIMD> p2 = mesh->Point (el[1]); + Point<3,SIMD> p3 = mesh->Point (el[2]); + + Vec<3,SIMD> vx = p1-p3; + Vec<3,SIMD> vy = p2-p3; + for (size_t ii = 0; ii < simd_npt; ii++) + { + simd_points[ii] = p3 + simd_pref[ii](0) * vx + simd_pref[ii](1) * vy; + for (size_t j = 0; j < 3; j++) + { + simd_dxdxis[ii](j,0) = vx(j); + simd_dxdxis[ii](j,1) = vy(j); + } + } + + Vec<3,SIMD> nv = Cross (vx, vy).Normalize(); + for (size_t ii = 0; ii < simd_npt; ii++) + simd_nvs[ii] = nv; + } + + + bool drawelem = false; + if (sol && sol->draw_surface) + { + // NgProfiler::StopTimer(timer1a); + // NgProfiler::StartTimer(timer1b); + drawelem = sol->solclass->GetMultiSurfValue (sei, -1, simd_npt, + &simd_pref[0](0).Data(), + &simd_points[0](0).Data(), + &simd_dxdxis[0](0).Data(), + &simd_values[0].Data()); + // NgProfiler::StopTimer(timer1b); + // NgProfiler::StartTimer(timer1c); + + for (size_t j = 0; j < sol->components; j++) + for (size_t i = 0; i < npt; i++) + mvalues[i*sol->components+j] = ((double*)&simd_values[j*simd_npt])[i]; + + if (usetexture == 2) + for (int ii = 0; ii < npt; ii++) + valuesc[ii] = ExtractValueComplex(sol, scalcomp, &mvalues[ii*sol->components]); + else + for (int ii = 0; ii < npt; ii++) + values[ii] = ExtractValue(sol, scalcomp, &mvalues[ii*sol->components]); + } + + for (size_t i = 0; i < npt; i++) + { + size_t ii = i/4; + size_t r = i%4; + for (int j = 0; j < 2; j++) + pref[i](j) = simd_pref[ii](j)[r]; + for (int j = 0; j < 3; j++) + points[i](j) = simd_points[ii](j)[r]; + for (int j = 0; j < 3; j++) + nvs[i](j) = simd_nvs[ii](j)[r]; + } + + if (deform) + for (int ii = 0; ii < npt; ii++) + points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1)); + + // NgProfiler::StopTimer(timer1c); + +#else bool curved = curv.IsSurfaceElementCurved(sei); for (int iy = 0, ii = 0; iy <= n; iy++) @@ -1492,6 +1677,8 @@ namespace netgen if (deform) for (int ii = 0; ii < npt; ii++) points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1)); +#endif + NgProfiler::StopTimer(timer1); int save_usetexture = usetexture; if (!drawelem) @@ -1500,108 +1687,62 @@ namespace netgen SetTextureMode (usetexture); } - for (int iy = 0, ii = 0; iy < n; iy++) + NgProfiler::StartTimer(timer2); + +#ifdef USE_BUFFERS + if (drawelem && usetexture == 1 && !logscale) { - glBegin (GL_TRIANGLE_STRIP); - for (int ix = 0; ix <= n-iy; ix++, ii++) - for (int k = 0; k < 2; k++) - { - if (ix+iy+k > n) continue; - int hi = (k == 0) ? ii : ii+n-iy+1; - - if (drawelem) - { - if (usetexture != 2) - SetOpenGlColor (values[hi]); - else - glTexCoord2f ( valuesc[hi].real(), valuesc[hi].imag() ); - } - else - glColor4fv (col_grey); - - glNormal3dv (nvs[hi]); - glVertex3dv (points[hi]); - } - glEnd(); + glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[0]); + glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(Point<3,double>), + npt*sizeof(Point<3,double>), &points[0][0]); + glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[1]); + glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(Vec<3,double>), + npt*sizeof(Vec<3,double>), &nvs[0][0]); + glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[2]); + glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(double), + npt*sizeof(double), &values[0]); + + for (size_t i = 0; i < ind_reftrig.Size(); i++) + glob_ind[i] = base_pi+ind_reftrig[i]; + + glBindBuffer (GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]); + glBufferSubData (GL_ELEMENT_ARRAY_BUFFER, surfel_vbo_size*sizeof(int), + ind_reftrig.Size()*sizeof(int), &glob_ind[0]); + surfel_vbo_size += ind_reftrig.Size(); + base_pi += npt; } + + else +#endif + for (int iy = 0, ii = 0; iy < n; iy++) + { + glBegin (GL_TRIANGLE_STRIP); + for (int ix = 0; ix <= n-iy; ix++, ii++) + for (int k = 0; k < 2; k++) + { + if (ix+iy+k > n) continue; + int hi = (k == 0) ? ii : ii+n-iy+1; + if (drawelem) + { + if (usetexture != 2) + SetOpenGlColor (values[hi]); + else + glTexCoord2f ( valuesc[hi].real(), valuesc[hi].imag() ); + } + else + glColor4fv (col_grey); + + glNormal3dv (nvs[hi]); + glVertex3dv (points[hi]); + } + glEnd(); + } + + NgProfiler::StopTimer(timer2); - /* - - GLuint vboId[3]; - glGenBuffersARB (3, &vboId[0]); -// cout << "vboId = " << vboId << endl; - - glBindBufferARB (GL_ARRAY_BUFFER_ARB, vboId[0]); - glBufferDataARB (GL_ARRAY_BUFFER_ARB, points.Size()*sizeof(Point<3>), - &points[0][0], GL_STATIC_DRAW_ARB); - - - // not so fast as old-fashened style - glEnableClientState(GL_VERTEX_ARRAY); - // glVertexPointer(3, GL_DOUBLE, 0, &points[0][0]); - glVertexPointer(3, GL_DOUBLE, 0, 0); //ARB - - glBindBufferARB (GL_ARRAY_BUFFER_ARB, vboId[1]); - glBufferDataARB (GL_ARRAY_BUFFER_ARB, nvs.Size()*sizeof(Point<3>), - &nvs[0][0], GL_STATIC_DRAW_ARB); - - glEnableClientState(GL_NORMAL_ARRAY); - // glNormalPointer(GL_DOUBLE, 0, &nvs[0][0]); - glNormalPointer(GL_DOUBLE, 0, 0); // ARB - - // if (drawelem && usetexture == 1) - { - glEnableClientState(GL_TEXTURE_COORD_ARRAY); - // glTexCoordPointer(1, GL_DOUBLE, 0, &values[0]); - - glBindBufferARB (GL_ARRAY_BUFFER_ARB, vboId[2]); - glBufferDataARB (GL_ARRAY_BUFFER_ARB, values.Size()*sizeof(double), - &values[0], GL_STATIC_DRAW_ARB); - glTexCoordPointer(1, GL_DOUBLE, 0, 0); - } - - Array gind; - - for (int iy = 0, ii = 0; iy < n; iy++,ii++) - { - for (int ix = 0; ix < n-iy; ix++, ii++) - { - int nv = (ix+iy+1 < n) ? 6 : 3; - - int ind[] = { ii, ii+1, ii+n-iy+1, - ii+n-iy+1, ii+1, ii+n-iy+2 }; - -// if (ix == 0 && iy == 0) -// for (int l = 0; l < 3; l++) -// { -// if (drawelem) -// { -// if (usetexture != 2) -// // SetOpenGlColor (values[ind[l]]); -// glTexCoord1f ( values[ind[l]] ); -// else -// glTexCoord2f ( valuesc[ind[l]].real(), valuesc[ind[l]].imag() ); -// } -// else -// glColor3fv (col_grey); -// } - - for (int j = 0; j < nv; j++) - gind.Append(ind[j]); - // glDrawElements(GL_TRIANGLES, nv, GL_UNSIGNED_INT, &ind[0]); - } - } - glDrawElements(GL_TRIANGLES, gind.Size(), GL_UNSIGNED_INT, &gind[0]); - - glDisableClientState(GL_VERTEX_ARRAY); - glDisableClientState(GL_NORMAL_ARRAY); - glDisableClientState(GL_TEXTURE_COORD_ARRAY); - - glDeleteBuffersARB (3, &vboId[0]); - */ - + if (!drawelem && (usetexture != save_usetexture)) { usetexture = save_usetexture; @@ -1609,9 +1750,32 @@ namespace netgen } } } + NgProfiler::StopTimer(timerloops); + + NgProfiler::StartTimer(timerbuffer); + + // glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]); + // glBufferData(GL_ELEMENT_ARRAY_BUFFER, glob_ind.Size()*sizeof(int), &glob_ind[0], GL_STATIC_DRAW); + // surfel_vbo_size = glob_ind.Size(); + + NgProfiler::StopTimer(timerbuffer); + + // glDrawElements(GL_TRIANGLES, surfel_vbo_size, GL_UNSIGNED_INT, 0); + + // glDrawElements(GL_TRIANGLES, glob_ind.Size(), GL_UNSIGNED_INT, &glob_ind[0]); + + // glDisableClientState(GL_VERTEX_ARRAY); + // glDisableClientState(GL_NORMAL_ARRAY); + // glDisableClientState(GL_TEXTURE_COORD_ARRAY); + + // glDeleteBuffers (1, &IndexVBOID); + // glDeleteBuffers (4, &vboId[0]); + + + NgProfiler::StartTimer(timerlist); glEndList (); - - + NgProfiler::StopTimer(timerlist); + #ifdef PARALLELGL glFinish(); if (id > 0) @@ -2492,12 +2656,14 @@ namespace netgen if(minv_local < minv) { lock_guard guard(min_mutex); - minv = minv_local; + if(minv_local < minv) + minv = minv_local; } if(maxv_local > maxv) { lock_guard guard(max_mutex); - maxv = maxv_local; + if(maxv_local > maxv) + maxv = maxv_local; } }); } diff --git a/libsrc/visualization/vssolution.hpp b/libsrc/visualization/vssolution.hpp index cfd31511..d8193d56 100644 --- a/libsrc/visualization/vssolution.hpp +++ b/libsrc/visualization/vssolution.hpp @@ -43,7 +43,12 @@ class DLL_HEADER VisualSceneSolution : public VisualScene Point<3> p; }; - + // #define USE_BUFFERS +#ifdef USE_BUFFERS + bool has_surfel_vbo = false; + GLuint surfel_vbo[4]; // + size_t surfel_vbo_size; +#endif int surfellist; int linelist; int element1dlist; diff --git a/ng/netgenpy.cpp b/ng/netgenpy.cpp index 9478dadf..32d9c694 100644 --- a/ng/netgenpy.cpp +++ b/ng/netgenpy.cpp @@ -1,7 +1,7 @@ // a wrapper to load netgen-dll into python #include -#include +#include <../general/ngpython.hpp> #ifdef WIN32 #define DLL_HEADER __declspec(dllimport) @@ -10,19 +10,26 @@ #endif -void DLL_HEADER ExportNetgenMeshing(); -void DLL_HEADER ExportMeshVis(); -void DLL_HEADER ExportCSG(); -void DLL_HEADER ExportCSGVis(); -void DLL_HEADER ExportGeom2d(); +void DLL_HEADER ExportNetgenMeshing(py::module &m); +void DLL_HEADER ExportMeshVis(py::module &m); +void DLL_HEADER ExportCSG(py::module &m); +void DLL_HEADER ExportCSGVis(py::module &m); +void DLL_HEADER ExportGeom2d(py::module &m); -BOOST_PYTHON_MODULE(libngpy) +PYBIND11_PLUGIN(libngpy) { - ExportCSG(); - ExportCSGVis(); - ExportNetgenMeshing(); - ExportMeshVis(); - ExportGeom2d(); + py::module ngpy("libngpy", "pybind netgen module"); + py::module csg = ngpy.def_submodule("_csg", "pybind csg module"); + ExportCSG(csg); + py::module csgvis = ngpy.def_submodule("csgvis", "pybind csgvis module"); + ExportCSGVis(csgvis); + py::module meshing = ngpy.def_submodule("_meshing", "pybind meshing module"); + ExportNetgenMeshing(meshing); + py::module meshvis = ngpy.def_submodule("meshvis", "pybind meshvis module"); + ExportMeshVis(meshvis); + py::module geom2d = ngpy.def_submodule("_geom2d", "pybind geom2d module"); + ExportGeom2d(geom2d); + return ngpy.ptr(); } // Force linking libnglib to libnetgenpy diff --git a/ng/ngvisual.tcl b/ng/ngvisual.tcl index 93a5618a..37cae1d2 100644 --- a/ng/ngvisual.tcl +++ b/ng/ngvisual.tcl @@ -897,8 +897,8 @@ proc fieldlinesdialog { } { set g $w.nb.settings - frame $g.linesettings -relief groove -borderwidth 3 - label $g.linesettings.title -text "\nLine Settings\n" + ttk::frame $g.linesettings -relief groove -borderwidth 3 + ttk::label $g.linesettings.title -text "\nLine Settings\n" # tixControl $g.linesettings.length -label "rel. Length: " -integer false \ # -variable visoptions.fieldlineslength -min 0.00001 -max 10000 -step 0.1 \ # -options { @@ -951,8 +951,8 @@ proc fieldlinesdialog { } { global visoptions.fieldlinestolerance - frame $g.odesettings -relief groove -borderwidth 3 - label $g.odesettings.title -text "\nODE Settings\n" + ttk::frame $g.odesettings -relief groove -borderwidth 3 + ttk::label $g.odesettings.title -text "\nODE Settings\n" # tixControl $g.odesettings.tol -label "rel. Tolerance: " -integer false \ # -variable visoptions.fieldlinestolerance -min 0.00001 -max 1 -step 0.01 \ # -options { diff --git a/ng/parallelfunc.cpp b/ng/parallelfunc.cpp index 2efb0259..d3e6d7d7 100644 --- a/ng/parallelfunc.cpp +++ b/ng/parallelfunc.cpp @@ -138,7 +138,7 @@ void ParallelRun() string redraw_cmd; // MyMPI_Recv (redraw_cmd, 0, MPI_TAG_VIS); redraw_cmd = MyMPI_RecvCmd(); - + // PrintMessage (1, "Redraw - ", redraw_cmd); static string displname; @@ -147,17 +147,17 @@ void ParallelRun() static Display * display = NULL; static GLXContext context; static XVisualInfo * visinfo = 0; - + // if (!display) if (redraw_cmd == "init") { MyMPI_Recv (displname, 0, MPI_TAG_VIS); MyMPI_Recv (curDrawable, 0, MPI_TAG_VIS); MyMPI_Recv (contextid, 0, MPI_TAG_VIS); - + display = XOpenDisplay (displname.c_str()); - + /* PrintMessage (3, "displ - name = ", displname); PrintMessage (3, "display = ", display, @@ -166,23 +166,27 @@ void ParallelRun() " , h = ", XDisplayHeight (display, 0)); */ - Window win; - int wx, wy; - unsigned int ww, wh, bw, depth; - // cout << "got drawable: " << curDrawable << endl; - - XGetGeometry(display, curDrawable, &win, - &wx, &wy, &ww, &wh, - &bw, &depth); - /* + Window win; + int wx, wy; + unsigned int ww, wh, bw, depth; + cout << "got drawable: " << curDrawable << ", contextid: " << contextid << endl; + + cout << "get geometriy..." << endl; + XGetGeometry(display, curDrawable, &win, + &wx, &wy, &ww, &wh, + &bw, &depth); + cout << "have!" << endl; + cout << "P" << id << ": window-props: x = " << wx << ", y = " << wy << ", w = " << ww << ", h = " << wh << ", depth = " << depth << endl; */ - + #define VISUAL #ifdef VISUAL - + +#ifdef VISINFO_OLD + //this does not seem to work anymore (but might still be with togl1.7?) // make a new GLXContext // first, generate a visual (copied from togl) @@ -221,51 +225,68 @@ void ParallelRun() attrib_list[attrib_count++] = GLX_DOUBLEBUFFER; attrib_list[attrib_count++] = None; - + visinfo = glXChooseVisual(display, 0, attrib_list); + cout << "have vis?" << endl; + if (visinfo) { /* found a GLX visual! */ // cout << "found VISINFO !!!" << endl; + cout << "found VISINFO !!!" << endl; /* - int hi = 0; - std::cout << "attribs = "; - while (attrib_list[hi] != None) + int hi = 0; + std::cout << "attribs = "; + while (attrib_list[hi] != None) std::cout << attrib_list[hi++] << " "; - std::cout << std::endl; + std::cout << std::endl; */ - + break; } } if (!visinfo) cerr << "no VISINFO found" << endl; +#else + //get all possible confs + int nconfs; + auto cptr = glXGetFBConfigs (display,0, &nconfs); + Array conf_ids(nconfs); + for(int k=0;kFBConfig->visual + unsigned int d_fbc_id; + glXQueryDrawable( display, curDrawable, GLX_FBCONFIG_ID, &d_fbc_id); + GLXFBConfig d_fbc; + for(int k=0;k