ForkMaster/netgen
1
0
Fork 0
mirror of https://github.com/NGSolve/netgen.git synced 2025-01-27 13:20:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into remove_tix2

Browse Source
This commit is contained in:
Gerhard Kitzler 2016-11-10 13:26:18 +01:00
commit fdd2fab64a
43 changed files with 1393 additions and 860 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitlab-ci.yml
View File

@ -136,7 +136,7 @@ build_ubuntu_1604:
paths: paths:
- build/ - build/
- src/ - src/
key: "netgen_win32" key: "netgen_win32_${CI_BUILD_REF_NAME}"
test_netgen_win64: test_netgen_win64:
<<: *win64 <<: *win64
@ -145,7 +145,7 @@ test_netgen_win64:
paths: paths:
- build/ - build/
- src/ - src/
key: "netgen_win64" key: "netgen_win64_${CI_BUILD_REF_NAME}"
# Linux # Linux
.template_test_linux: &test_linux .template_test_linux: &test_linux

3
.gitmodules vendored Normal file
View File

@ -0,0 +1,3 @@
[submodule "external_dependencies/pybind11"]
path = external_dependencies/pybind11
url = https://github.com/pybind/pybind11.git

54
CMakeLists.txt
View File

@ -124,7 +124,7 @@ endmacro()
if(WIN32) if(WIN32)
get_WIN32_WINNT(ver) 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) set(CMAKE_MFC_FLAG 0)
# add_definitions(-DNGINTERFACE_EXPORTS) # add_definitions(-DNGINTERFACE_EXPORTS)
# add_definitions(-DNGLIB_EXPORTS) # add_definitions(-DNGLIB_EXPORTS)
@ -135,7 +135,7 @@ if(WIN32)
else(WIN32) else(WIN32)
# build shared libraries # build shared libraries
set(NG_LIB_TYPE SHARED) 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) endif(WIN32)
if(APPLE) if(APPLE)
# set(MACOSX_BUNDLE ON) # set(MACOSX_BUNDLE ON)
@ -185,6 +185,17 @@ endif (USE_GUI)
####################################################################### #######################################################################
if (USE_PYTHON) 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(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" "${PROJECT_SOURCE_DIR}/cmake_modules/python")
set(PYTHON_VERSION "3" CACHE STRING "") set(PYTHON_VERSION "3" CACHE STRING "")
set(Python_ADDITIONAL_VERSIONS 3.5) set(Python_ADDITIONAL_VERSIONS 3.5)
@ -198,36 +209,10 @@ if (USE_PYTHON)
find_package(PythonInterp ${PYTHON_VERSION} EXACT REQUIRED) find_package(PythonInterp ${PYTHON_VERSION} EXACT REQUIRED)
find_package(PythonLibs ${PYTHON_VERSION} EXACT REQUIRED) find_package(PythonLibs ${PYTHON_VERSION} EXACT REQUIRED)
endif( PYTHON_VERSION STREQUAL "3") 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) 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(${PYTHON_INCLUDE_DIRS})
include_directories(${Boost_INCLUDE_DIRS}) set(PYTHON_LIBS "${PYTHON_LIBRARIES}")
set(PYTHON_LIBS "${PYTHON_LIBRARIES};${Boost_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) 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}") message("python install path: ${PYTHON_PACKAGES_INSTALL_DIR}")
endif (USE_PYTHON) endif (USE_PYTHON)
@ -288,7 +273,7 @@ if(INSTALL_DEPENDENCIES)
# find_library(LIBSTDCXX NAMES stdc++.6) # find_library(LIBSTDCXX NAMES stdc++.6)
# find_library(LIBGCCS NAMES gcc_s.1) # 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}" ${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) get_filename_component(MY_LIB_DIR ${TK_LIBRARY} DIRECTORY)
if(APPLE) if(APPLE)
if(TIX_LIBRARY) if(TIX_LIBRARY)
@ -355,8 +340,7 @@ if(UNIX)
string(SUBSTRING ${temp} 17 -1 UBUNTU_VERSION) string(SUBSTRING ${temp} 17 -1 UBUNTU_VERSION)
message("ubuntu version: ${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, libtk8.5, libtcl8.5, tix, libxmu6")
set(CPACK_DEBIAN_PACKAGE_DEPENDS "python3, ${BOOST_PACKAGE}, libtk8.5, libtcl8.5, tix, libxmu6")
execute_process(COMMAND dpkg --print-architecture OUTPUT_VARIABLE CPACK_DEBIAN_PACKAGE_ARCHITECTURE OUTPUT_STRIP_TRAILING_WHITESPACE) execute_process(COMMAND dpkg --print-architecture OUTPUT_VARIABLE CPACK_DEBIAN_PACKAGE_ARCHITECTURE OUTPUT_STRIP_TRAILING_WHITESPACE)
set(CPACK_DEBIAN_PACKAGE_MAINTAINER "Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>") set(CPACK_DEBIAN_PACKAGE_MAINTAINER "Matthias Hochsteger <matthias.hochsteger@tuwien.ac.at>")
if(USE_MPI) if(USE_MPI)
@ -364,11 +348,7 @@ if(UNIX)
set(CPACK_PACKAGE_NAME "${CPACK_PACKAGE_NAME}_mpi") set(CPACK_PACKAGE_NAME "${CPACK_PACKAGE_NAME}_mpi")
endif(USE_MPI) endif(USE_MPI)
if(USE_OCC) if(USE_OCC)
if("${UBUNTU_VERSION}" STREQUAL "xenial") set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf-dev")
set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf10")
else()
set(CPACK_DEBIAN_PACKAGE_DEPENDS "${CPACK_DEBIAN_PACKAGE_DEPENDS}, liboce-ocaf8")
endif()
endif(USE_OCC) endif(USE_OCC)
set(CPACK_DEBIAN_PACKAGE_SECTION Science) set(CPACK_DEBIAN_PACKAGE_SECTION Science)
set(CPACK_DEBIAN_PACKAGE_NAME ${CPACK_PACKAGE_NAME}) set(CPACK_DEBIAN_PACKAGE_NAME ${CPACK_PACKAGE_NAME})

1
external_dependencies/pybind11 Submodule

@ -0,0 +1 @@
Subproject commit 2b92a49115657712c7dd924248df2286e6143b8d

2
libsrc/csg/csgeom.cpp
View File

@ -1494,7 +1494,7 @@ namespace netgen
// virtual VisualScene * GetVisualScene (const NetgenGeometry * geom) const; // 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 NetgenGeometry * CSGeometryRegister :: Load (string filename) const
{ {

10
libsrc/csg/csgparser.cpp
View File

@ -822,11 +822,17 @@ namespace netgen
/* /*
Main parsing function for CSG geometry Main parsing function for CSG geometry
*/ */
CSGeometry * ParseCSG (istream & istr) CSGeometry * ParseCSG (istream & istr, CSGeometry * instance=nullptr)
{ {
CSGScanner scan(istr); CSGScanner scan(istr);
geom = new CSGeometry; if (instance)
{
new (instance) CSGeometry;
geom = instance;
}
else
geom = new CSGeometry;
scan.ReadNext(); scan.ReadNext();
if (scan.GetToken() != TOK_RECO) // keyword 'algebraic3d' if (scan.GetToken() != TOK_RECO) // keyword 'algebraic3d'

2
libsrc/csg/genmesh.cpp
View File

@ -687,6 +687,8 @@ namespace netgen
mparam.grading); mparam.grading);
mesh -> LoadLocalMeshSize (mparam.meshsizefilename); mesh -> LoadLocalMeshSize (mparam.meshsizefilename);
for (auto mspnt : mparam.meshsize_points)
mesh -> RestrictLocalH (mspnt.pnt, mspnt.h);
} }
spoints.SetSize(0); spoints.SetSize(0);

12
libsrc/csg/identify.cpp
View File

@ -1079,6 +1079,7 @@ void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh)
for (int i = 1; i <= mesh.GetNFD(); i++) for (int i = 1; i <= mesh.GetNFD(); i++)
{ {
auto & fdi = mesh.GetFaceDescriptor(i);
int surfi = mesh.GetFaceDescriptor(i).SurfNr(); int surfi = mesh.GetFaceDescriptor(i).SurfNr();
if (s1rep != surfi) continue; if (s1rep != surfi) continue;
@ -1090,12 +1091,21 @@ void CloseSurfaceIdentification :: IdentifyFaces (class Mesh & mesh)
for (int j = 1; j <= mesh.GetNFD(); j++) for (int j = 1; j <= mesh.GetNFD(); j++)
{ {
auto & fdj = mesh.GetFaceDescriptor(j);
int surfj = mesh.GetFaceDescriptor(j).SurfNr(); int surfj = mesh.GetFaceDescriptor(j).SurfNr();
if (surfi == surfj) continue; if (surfi == surfj) continue;
if (s2rep != 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; int idok = 1;
for (side = 1; side <= 2 && idok; side++) for (side = 1; side <= 2 && idok; side++)

248
libsrc/csg/python_csg.cpp
View File

@ -11,7 +11,6 @@ namespace netgen
extern shared_ptr<NetgenGeometry> ng_geometry; extern shared_ptr<NetgenGeometry> ng_geometry;
} }
inline void NOOP_Deleter(void *) { ; }
// a shadow solid tree using shared pointers. // a shadow solid tree using shared pointers.
@ -160,60 +159,103 @@ inline ostream & operator<< (ostream & ost, const SPSolid & sol)
namespace netgen 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"); py::class_<NGDummyArgument>(m, "NGDummyArgument")
.def("__bool__", []( NGDummyArgument &self ) { return false; } )
bp::class_<Point<2>> ("Point2d", bp::init<double,double>()) ;
py::class_<Point<2>> (m, "Point2d")
.def(py::init<double,double>())
.def ("__str__", &ToString<Point<2>>) .def ("__str__", &ToString<Point<2>>)
.def(bp::self-bp::self) .def(py::self-py::self)
.def(bp::self+Vec<2>()) .def(py::self+Vec<2>())
.def(bp::self-Vec<2>()) .def(py::self-Vec<2>())
; ;
bp::class_<Point<3>> ("Point3d", bp::init<double,double,double>()) py::class_<Point<3>> (m, "Point3d")
.def(py::init<double,double,double>())
.def ("__str__", &ToString<Point<3>>) .def ("__str__", &ToString<Point<3>>)
.def(bp::self-bp::self) .def(py::self-py::self)
.def(bp::self+Vec<3>()) .def(py::self+Vec<3>())
.def(bp::self-Vec<3>()) .def(py::self-Vec<3>())
; ;
bp::def ("Pnt", FunctionPointer m.def ("Pnt", FunctionPointer
([](double x, double y, double z) { return Point<3>(x,y,z); })); ([](double x, double y, double z) { return global_trafo(Point<3>(x,y,z)); }));
bp::def ("Pnt", FunctionPointer m.def ("Pnt", FunctionPointer
([](double x, double y) { return Point<2>(x,y); })); ([](double x, double y) { return Point<2>(x,y); }));
bp::class_<Vec<2>> ("Vec2d", bp::init<double,double>())
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_<Vec<2>> (m, "Vec2d")
.def(py::init<double,double>())
.def ("__str__", &ToString<Vec<3>>) .def ("__str__", &ToString<Vec<3>>)
.def(bp::self+bp::self) .def(py::self+py::self)
.def(bp::self-bp::self) .def(py::self-py::self)
.def(-bp::self) .def(-py::self)
.def(double()*bp::self) .def(double()*py::self)
.def("Norm", &Vec<2>::Length) .def("Norm", &Vec<2>::Length)
; ;
bp::class_<Vec<3>> ("Vec3d", bp::init<double,double,double>()) py::class_<Vec<3>> (m, "Vec3d")
.def(py::init<double,double,double>())
.def ("__str__", &ToString<Vec<3>>) .def ("__str__", &ToString<Vec<3>>)
.def(bp::self+bp::self) .def(py::self+py::self)
.def(bp::self-bp::self) .def(py::self-py::self)
.def(-bp::self) .def(-py::self)
.def(double()*bp::self) .def(double()*py::self)
.def("Norm", &Vec<3>::Length) .def("Norm", &Vec<3>::Length)
; ;
bp::def ("Vec", FunctionPointer m.def ("Vec", FunctionPointer
([] (double x, double y, double z) { return Vec<3>(x,y,z); })); ([] (double x, double y, double z) { return global_trafo(Vec<3>(x,y,z)); }));
bp::def ("Vec", FunctionPointer m.def ("Vec", FunctionPointer
([] (double x, double y) { return Vec<2>(x,y); })); ([] (double x, double y) { return Vec<2>(x,y); }));
#if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100) py::class_<SplineGeometry<2>> (m, "SplineCurve2d")
bp::register_ptr_to_python<shared_ptr<SPSolid>>(); .def(py::init<>())
#endif .def ("AddPoint", FunctionPointer
bp::class_<SPSolid, shared_ptr<SPSolid>, boost::noncopyable> ("Solid", bp::no_init) ([] (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_<SPSolid, shared_ptr<SPSolid>> (m, "Solid")
.def ("__str__", &ToString<SPSolid>) .def ("__str__", &ToString<SPSolid>)
.def ("__add__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::UNION, self, other); })) .def ("__add__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::UNION, self, other); }))
.def ("__mul__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::SECTION, self, other); })) .def ("__mul__", FunctionPointer( [] ( shared_ptr<SPSolid> self, shared_ptr<SPSolid> other) { return make_shared<SPSolid> (SPSolid::SECTION, self, other); }))
@ -231,100 +273,105 @@ DLL_HEADER void ExportCSG()
.def ("mat", FunctionPointer([](shared_ptr<SPSolid> & self, string mat) -> shared_ptr<SPSolid> .def ("mat", FunctionPointer([](shared_ptr<SPSolid> & self, string mat) -> shared_ptr<SPSolid>
{ self->SetMaterial(mat); return self; })) { self->SetMaterial(mat); return self; }))
.def ("mat", &SPSolid::GetMaterial) .def ("mat", &SPSolid::GetMaterial)
.def("col", FunctionPointer([](shared_ptr<SPSolid> & self, bp::list rgb) -> shared_ptr<SPSolid> .def("col", FunctionPointer([](shared_ptr<SPSolid> & self, py::list rgb) -> shared_ptr<SPSolid>
{ {
bp::extract<double> red(rgb[0]); py::extract<double> red(rgb[0]);
bp::extract<double> green(rgb[1]); py::extract<double> green(rgb[1]);
bp::extract<double> blue(rgb[2]); py::extract<double> blue(rgb[2]);
self->SetColor(red(),green(),blue()); self->SetColor(red(),green(),blue());
return self; return self;
})) }))
.def("transp", FunctionPointer([](shared_ptr<SPSolid> & self)->shared_ptr < SPSolid > { self->SetTransparent(); return self; })) .def("transp", FunctionPointer([](shared_ptr<SPSolid> & 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); Sphere * sp = new Sphere (c, r);
Solid * sol = new Solid (sp); Solid * sol = new Solid (sp);
return make_shared<SPSolid> (sol); return make_shared<SPSolid> (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); Plane * sp = new Plane (p,n);
Solid * sol = new Solid (sp); Solid * sol = new Solid (sp);
return make_shared<SPSolid> (sol); return make_shared<SPSolid> (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); Cone * cyl = new Cone (a, b, ra, rb);
Solid * sol = new Solid (cyl); Solid * sol = new Solid (cyl);
return make_shared<SPSolid> (sol); return make_shared<SPSolid> (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); Cylinder * cyl = new Cylinder (a, b, r);
Solid * sol = new Solid (cyl); Solid * sol = new Solid (cyl);
return make_shared<SPSolid> (sol); return make_shared<SPSolid> (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); OrthoBrick * brick = new OrthoBrick (p1,p2);
Solid * sol = new Solid (brick); Solid * sol = new Solid (brick);
return make_shared<SPSolid> (sol); return make_shared<SPSolid> (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<SPSolid> (sol);
}));
bp::def ("Or", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2) m.def ("Or", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
{ {
return make_shared<SPSolid> (SPSolid::UNION, s1, s2); return make_shared<SPSolid> (SPSolid::UNION, s1, s2);
})); }));
bp::def ("And", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2) m.def ("And", FunctionPointer([](shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2)
{ {
return make_shared<SPSolid> (SPSolid::SECTION, s1, s2); return make_shared<SPSolid> (SPSolid::SECTION, s1, s2);
})); }));
bp::class_<CSGeometry,shared_ptr<CSGeometry>, boost::noncopyable> ("CSGeometry") py::class_<CSGeometry,shared_ptr<CSGeometry>> (m, "CSGeometry")
.def("__init__", bp::make_constructor (FunctionPointer .def(py::init<>())
([](const string & filename) .def("__init__",
[](CSGeometry *instance, const string & filename)
{ {
cout << "load geometry"; cout << "load geometry";
ifstream ist(filename); ifstream ist(filename);
shared_ptr<CSGeometry> geom(ParseCSG(ist)); ParseCSG(ist, instance);
geom -> FindIdenticSurfaces(1e-8 * geom->MaxSize()); instance -> FindIdenticSurfaces(1e-8 * instance->MaxSize());
return geom; })
}))) .def("__init__",
[](CSGeometry *instance, const py::list & solidlist)
.def("__init__", bp::make_constructor (FunctionPointer
([](const bp::list & solidlist)
{ {
cout << "csg from list"; cout << "csg from list";
auto geom = make_shared<CSGeometry>(); new (instance) CSGeometry();
for (int i = 0; i < len(solidlist); i++) for (int i = 0; i < len(solidlist); i++)
{ {
bp::object obj = solidlist[i]; py::object obj = solidlist[i];
cout << "obj " << i << endl; cout << "obj " << i << endl;
bp::extract<shared_ptr<SPSolid>> solid(solidlist[i]); py::extract<shared_ptr<SPSolid>> solid(solidlist[i]);
if(solid.check()) if(solid.check())
{ {
cout << "its a solid" << endl; cout << "its a solid" << endl;
solid()->AddSurfaces (*geom); solid()->AddSurfaces (*instance);
solid()->GiveUpOwner(); solid()->GiveUpOwner();
int tlonr = geom->SetTopLevelObject (solid()->GetSolid()); int tlonr = instance->SetTopLevelObject (solid()->GetSolid());
geom->GetTopLevelObject(tlonr) -> SetMaterial(solid()->GetMaterial()); instance->GetTopLevelObject(tlonr) -> SetMaterial(solid()->GetMaterial());
} }
} }
geom -> FindIdenticSurfaces(1e-8 * geom->MaxSize()); instance -> FindIdenticSurfaces(1e-8 * instance->MaxSize());
return geom; })
})))
.def("Save", FunctionPointer([] (CSGeometry & self, string filename) .def("Save", FunctionPointer([] (CSGeometry & self, string filename)
{ {
cout << "save geometry to file " << filename << endl; cout << "save geometry to file " << filename << endl;
self.Save (filename); self.Save (filename);
})) }))
.def("Add", FunctionPointer .def("Add",
([] (CSGeometry & self, shared_ptr<SPSolid> solid, bp::list bcmod) [] (CSGeometry & self, shared_ptr<SPSolid> solid, py::list bcmod)
{ {
solid->AddSurfaces (self); solid->AddSurfaces (self);
solid->GiveUpOwner(); solid->GiveUpOwner();
@ -334,15 +381,15 @@ DLL_HEADER void ExportCSG()
self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent()); self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
// bcmod is list of tuples ( solid, bcnr ) // 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<bp::tuple> (bcmod[i]) (); py::tuple tup = py::extract<py::tuple> (bcmod[i]) ();
auto mod_solid = bp::extract<shared_ptr<SPSolid>> (tup[0]) (); auto mod_solid = py::extract<shared_ptr<SPSolid>> (tup[0]) ();
int mod_nr = -1; int mod_nr = -1;
string * bcname = nullptr; string * bcname = nullptr;
bp::object val = tup[1]; py::object val = tup[1];
if (bp::extract<int>(val).check()) mod_nr = bp::extract<int> (val)(); if (py::extract<int>(val).check()) mod_nr = py::extract<int> (val)();
if (bp::extract<string>(val).check()) bcname = new string ( bp::extract<string> (val)()); if (py::extract<string>(val).check()) bcname = new string ( py::extract<string> (val)());
Array<int> si; Array<int> si;
mod_solid -> GetSolid() -> GetSurfaceIndices (si); mod_solid -> GetSolid() -> GetSurfaceIndices (si);
@ -360,9 +407,8 @@ DLL_HEADER void ExportCSG()
delete bcname; delete bcname;
} }
return tlonr; return tlonr;
},
}), py::arg("solid"), py::arg("bcmod")=py::list()
(bp::arg("self"), bp::arg("solid"), bp::arg("bcmod")=bp::list())
) )
.def("AddSurface", FunctionPointer .def("AddSurface", FunctionPointer
@ -379,13 +425,13 @@ DLL_HEADER void ExportCSG()
self.GetTopLevelObject(tlonr) -> SetRGB(solid->GetRed(),solid->GetGreen(),solid->GetBlue()); self.GetTopLevelObject(tlonr) -> SetRGB(solid->GetRed(),solid->GetGreen(),solid->GetBlue());
self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent()); self.GetTopLevelObject(tlonr)->SetTransparent(solid->IsTransparent());
}), }),
(bp::arg("self"), bp::arg("surface"), bp::arg("solid")) py::arg("surface"), py::arg("solid")
) )
.def("CloseSurfaces", FunctionPointer .def("CloseSurfaces", FunctionPointer
([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, bp::list aslices ) ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, py::list aslices )
{ {
Array<int> si1, si2; Array<int> si1, si2;
s1->GetSolid()->GetSurfaceIndices (si1); s1->GetSolid()->GetSurfaceIndices (si1);
@ -397,15 +443,15 @@ DLL_HEADER void ExportCSG()
try try
{ {
int n = bp::len(aslices); int n = py::len(aslices);
Array<double> slices(n); Array<double> slices(n);
for(int i=0; i<n; i++) for(int i=0; i<n; i++)
{ {
slices[i]= bp::extract<double>(aslices[i])(); slices[i]= py::extract<double>(aslices[i])();
} }
flags.SetFlag("slices", slices); flags.SetFlag("slices", slices);
} }
catch( bp::error_already_set const & ) { catch( py::error_already_set const & ) {
cout << "caught python error:" << endl; cout << "caught python error:" << endl;
PyErr_Print(); PyErr_Print();
} }
@ -418,7 +464,7 @@ DLL_HEADER void ExportCSG()
domain, domain,
flags)); 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 .def("CloseSurfaces", FunctionPointer
([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, int reflevels) ([] (CSGeometry & self, shared_ptr<SPSolid> s1, shared_ptr<SPSolid> s2, int reflevels)
@ -438,7 +484,7 @@ DLL_HEADER void ExportCSG()
domain, domain,
flags)); 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 .def("PeriodicSurfaces", FunctionPointer
@ -453,7 +499,7 @@ DLL_HEADER void ExportCSG()
(self.GetNIdentifications()+1, self, (self.GetNIdentifications()+1, self,
self.GetSurface (si1[0]), self.GetSurface (si2[0]))); 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 .def("GetTransparent", FunctionPointer
@ -461,14 +507,14 @@ DLL_HEADER void ExportCSG()
{ {
return self.GetTopLevelObject(tlonr)->GetTransparent(); return self.GetTopLevelObject(tlonr)->GetTransparent();
}), }),
(bp::arg("self"), bp::arg("tlonr")) py::arg("tlonr")
) )
.def("SetTransparent", FunctionPointer .def("SetTransparent", FunctionPointer
([] (CSGeometry & self, int tlonr, bool transparent) ([] (CSGeometry & self, int tlonr, bool transparent)
{ {
self.GetTopLevelObject(tlonr)->SetTransparent(transparent); self.GetTopLevelObject(tlonr)->SetTransparent(transparent);
}), }),
(bp::arg("self"), bp::arg("tlonr"), bp::arg("transparent")) py::arg("tlonr"), py::arg("transparent")
) )
.def("GetVisible", FunctionPointer .def("GetVisible", FunctionPointer
@ -476,26 +522,34 @@ DLL_HEADER void ExportCSG()
{ {
return self.GetTopLevelObject(tlonr)->GetVisible(); return self.GetTopLevelObject(tlonr)->GetVisible();
}), }),
(bp::arg("self"), bp::arg("tlonr")) py::arg("tlonr")
) )
.def("SetVisible", FunctionPointer .def("SetVisible", FunctionPointer
([] (CSGeometry & self, int tlonr, bool visible) ([] (CSGeometry & self, int tlonr, bool visible)
{ {
self.GetTopLevelObject(tlonr)->SetVisible(visible); self.GetTopLevelObject(tlonr)->SetVisible(visible);
}), }),
(bp::arg("self"), bp::arg("tlonr"), bp::arg("visible")) py::arg("tlonr"), py::arg("visible")
) )
.def("SetBoundingBox", FunctionPointer .def("SetBoundingBox", FunctionPointer
([] (CSGeometry & self, Point<3> pmin, Point<3> pmax) ([] (CSGeometry & self, Point<3> pmin, Point<3> pmax)
{ {
self.SetBoundingBox(Box<3> (pmin, 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<CSGeometry> 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<CSGeometry> geo, MeshingParameters & param) ([](shared_ptr<CSGeometry> geo, MeshingParameters & param)
{ {
auto dummy = make_shared<Mesh>(); auto dummy = make_shared<Mesh>();
@ -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) ([](const Mesh & self, const string & filename, const CSGeometry & geom)
{ {
ostream * outfile; ostream * outfile;
@ -533,7 +587,7 @@ DLL_HEADER void ExportCSG()
bp::def("ZRefinement", FunctionPointer m.def("ZRefinement", FunctionPointer
([](Mesh & mesh, CSGeometry & geom) ([](Mesh & mesh, CSGeometry & geom)
{ {
ZRefinementOptions opt; ZRefinementOptions opt;
@ -543,14 +597,10 @@ DLL_HEADER void ExportCSG()
; ;
} }
PYBIND11_PLUGIN(libcsg) {
py::module m("csg", "pybind csg");
ExportCSG(m);
return m.ptr();
BOOST_PYTHON_MODULE(libcsg) {
ExportCSG();
} }
#endif #endif

73
libsrc/csg/revolution.cpp
View File

@ -156,25 +156,42 @@ namespace netgen
{ {
if(spline_coefficient.Size() == 0) if(spline_coefficient.Size() == 0)
spline->GetCoeff(spline_coefficient); spline->GetCoeff(spline_coefficient);
if(spline_coefficient_shifted.Size() == 0)
spline->GetCoeff(spline_coefficient_shifted, spline->StartPI());
Point<2> p; Point<2> p;
CalcProj(point,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(2)*p(0)*p(1) + spline_coefficient(3)*p(0)
+ spline_coefficient(4)*p(1) + spline_coefficient(5); + 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 void RevolutionFace :: CalcGradient (const Point<3> & point, Vec<3> & grad) const
{ {
if(spline_coefficient.Size() == 0) if(spline_coefficient.Size() == 0)
spline->GetCoeff(spline_coefficient); spline->GetCoeff(spline_coefficient);
if(spline_coefficient_shifted.Size() == 0)
spline->GetCoeff(spline_coefficient_shifted, spline->StartPI());
Vec<3> point_minus_p0 = point-p0; Vec<3> point_minus_p0 = point-p0;
Point<2> p; Point<2> p;
CalcProj0(point_minus_p0,p); CalcProj0(point_minus_p0,p);
/*
const double dFdxbar = 2.*spline_coefficient(0)*p(0) + spline_coefficient(2)*p(1) + spline_coefficient(3); const double dFdxbar = 2.*spline_coefficient(0)*p(0) + spline_coefficient(2)*p(1) + spline_coefficient(3);
if(fabs(p(1)) > 1e-10) if(fabs(p(1)) > 1e-10)
@ -193,6 +210,27 @@ namespace netgen
grad(2) = dFdxbar*v_axis(2); grad(2) = dFdxbar*v_axis(2);
//(*testout) << "grad2("<<point<<") = " << grad << endl; //(*testout) << "grad2("<<point<<") = " << grad << endl;
} }
*/
Vec<2> 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("<<point<<") = " << grad << endl;
}
else
{
grad(0) = dFdxbar*v_axis(0);
grad(1) = dFdxbar*v_axis(1);
grad(2) = dFdxbar*v_axis(2);
//(*testout) << "grad2("<<point<<") = " << grad << endl;
}
} }
@ -551,21 +589,35 @@ namespace netgen
} }
*/ */
INSOLID_TYPE RevolutionFace :: PointInFace (const Point<3> & p, const double eps) const /* INSOLID_TYPE */ bool RevolutionFace :: PointInFace (const Point<3> & p, const double eps) const
{ {
Point<2> p2d; Point<2> p2d;
CalcProj(p,p2d); CalcProj(p,p2d);
if (!spline -> InConvexHull(p2d, eps)) return false;
/*
double val = spline_coefficient(0)*p2d(0)*p2d(0) + spline_coefficient(1)*p2d(1)*p2d(1) + spline_coefficient(2)*p2d(0)*p2d(1) + double val = spline_coefficient(0)*p2d(0)*p2d(0) + spline_coefficient(1)*p2d(1)*p2d(1) + spline_coefficient(2)*p2d(0)*p2d(1) +
spline_coefficient(3)*p2d(0) + spline_coefficient(4)*p2d(1) + spline_coefficient(5); spline_coefficient(3)*p2d(0) + spline_coefficient(4)*p2d(1) + spline_coefficient(5);
*/
Vec<2> pr = p2d - spline->StartPI();
double val = 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);
return (fabs(val) < eps);
/*
if(val > eps) if(val > eps)
return IS_OUTSIDE; return IS_OUTSIDE;
if(val < -eps) if(val < -eps)
return IS_INSIDE; return IS_INSIDE;
return DOES_INTERSECT; return DOES_INTERSECT;
*/
} }
@ -748,6 +800,15 @@ namespace netgen
return IS_INSIDE; return IS_INSIDE;
} }
void Revolution :: GetTangentialSurfaceIndices (const Point<3> & p,
Array<int> & surfind, double eps) const
{
for (int j = 0; j < faces.Size(); j++)
if (faces[j] -> PointInFace(p, eps))
if (!surfind.Contains (GetSurfaceId(j)))
surfind.Append (GetSurfaceId(j));
}
INSOLID_TYPE Revolution :: VecInSolid (const Point<3> & p, INSOLID_TYPE Revolution :: VecInSolid (const Point<3> & p,
const Vec<3> & v, const Vec<3> & v,
double eps) const double eps) const
@ -763,7 +824,7 @@ namespace netgen
Array<int> intersecting_faces; Array<int> intersecting_faces;
for(int i=0; i<faces.Size(); i++) for(int i=0; i<faces.Size(); i++)
if(faces[i]->PointInFace(p,eps) == DOES_INTERSECT) if(faces[i]->PointInFace(p,eps)) // == DOES_INTERSECT)
intersecting_faces.Append(i); intersecting_faces.Append(i);
Vec<3> hv; Vec<3> hv;

10
libsrc/csg/revolution.hpp
View File

@ -17,8 +17,10 @@ namespace netgen
Vec<3> v_axis; Vec<3> v_axis;
int id; int id;
// coefficient for implicizt polynomial
mutable Vector spline_coefficient; mutable Vector spline_coefficient;
mutable Vector spline_coefficient_shifted;
Array < Vec<2>* > checklines_vec; Array < Vec<2>* > checklines_vec;
@ -74,7 +76,7 @@ namespace netgen
const SplineSeg<2> & GetSpline(void) const {return *spline;} const SplineSeg<2> & GetSpline(void) const {return *spline;}
INSOLID_TYPE PointInFace (const Point<3> & p, const double eps) const; /* INSOLID_TYPE */ bool PointInFace (const Point<3> & p, const double eps) const;
void GetRawData(Array<double> & data) const; void GetRawData(Array<double> & data) const;
@ -125,6 +127,10 @@ namespace netgen
virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const; virtual INSOLID_TYPE BoxInSolid (const BoxSphere<3> & box) const;
virtual INSOLID_TYPE PointInSolid (const Point<3> & p, virtual INSOLID_TYPE PointInSolid (const Point<3> & p,
double eps) const; double eps) const;
virtual void GetTangentialSurfaceIndices (const Point<3> & p,
Array<int> & surfind, double eps) const;
virtual INSOLID_TYPE VecInSolid (const Point<3> & p, virtual INSOLID_TYPE VecInSolid (const Point<3> & p,
const Vec<3> & v, const Vec<3> & v,
double eps) const; double eps) const;

32
libsrc/csg/vscsg.cpp
View File

@ -488,25 +488,21 @@ namespace netgen
#ifdef NG_PYTHON #ifdef NG_PYTHON
#include <boost/python.hpp>
#include <../general/ngpython.hpp> #include <../general/ngpython.hpp>
namespace bp = boost::python;
DLL_HEADER void ExportCSGVis() DLL_HEADER void ExportCSGVis(py::module &m)
{ {
using namespace netgen; using namespace netgen;
ModuleScope module("csgvis"); py::class_<VisualSceneGeometry, shared_ptr<VisualSceneGeometry>>
(m, "VisualSceneGeometry")
bp::class_<VisualSceneGeometry, shared_ptr<VisualSceneGeometry>>
("VisualSceneGeometry", bp::no_init)
.def("Draw", &VisualSceneGeometry::DrawScene) .def("Draw", &VisualSceneGeometry::DrawScene)
; ;
bp::def("SetBackGroundColor", &VisualSceneGeometry::SetBackGroundColor); m.def("SetBackGroundColor", &VisualSceneGeometry::SetBackGroundColor);
bp::def("VS", FunctionPointer m.def("VS",
([](CSGeometry & geom) [](CSGeometry & geom)
{ {
geom.FindIdenticSurfaces(1e-6); geom.FindIdenticSurfaces(1e-6);
geom.CalcTriangleApproximation(0.01, 20); geom.CalcTriangleApproximation(0.01, 20);
@ -514,17 +510,19 @@ DLL_HEADER void ExportCSGVis()
vs->SetGeometry(&geom); vs->SetGeometry(&geom);
return vs; return vs;
})); });
bp::def("MouseMove", FunctionPointer m.def("MouseMove",
([](VisualSceneGeometry &vsgeom, int oldx, int oldy, int newx, int newy, char mode) [](VisualSceneGeometry &vsgeom, int oldx, int oldy, int newx, int newy, char mode)
{ {
vsgeom.MouseMove(oldx, oldy, newx, newy, mode); vsgeom.MouseMove(oldx, oldy, newx, newy, mode);
})); });
} }
BOOST_PYTHON_MODULE(libcsgvis)
{ PYBIND11_PLUGIN(libcsgvis) {
ExportCSGVis(); py::module m("csg", "pybind csg");
ExportCSGVis(m);
return m.ptr();
} }
#endif #endif

2
libsrc/general/hashtabl.hpp
View File

@ -779,7 +779,7 @@ protected:
mask = size-1; mask = size-1;
// cout << "mask = " << mask << endl; // cout << "mask = " << mask << endl;
invalid = -1; invalid = -1;
for (int i = 0; i < size; i++) for (size_t i = 0; i < size; i++)
hash[i].I1() = invalid; hash[i].I1() = invalid;
} }

8
libsrc/general/mpi_interface.cpp
View File

@ -16,17 +16,17 @@ namespace netgen
void MyMPI_SendCmd (const char * cmd) void MyMPI_SendCmd (const char * cmd)
{ {
char buf[100]; char buf[10000];
strcpy (buf, cmd); strcpy (buf, cmd);
// MPI_Bcast (&buf, 100, MPI_CHAR, 0, MPI_COMM_WORLD); // MPI_Bcast (&buf, 100, MPI_CHAR, 0, MPI_COMM_WORLD);
for (int dest = 1; dest < ntasks; dest++) 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 () string MyMPI_RecvCmd ()
{ {
char buf[100]; char buf[10000];
// MPI_Bcast (&buf, 100, MPI_CHAR, 0, MPI_COMM_WORLD); // MPI_Bcast (&buf, 100, MPI_CHAR, 0, MPI_COMM_WORLD);
// VT_OFF(); // VT_OFF();
@ -44,7 +44,7 @@ namespace netgen
while (!flag); while (!flag);
// VT_ON(); // 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); return string(buf);
} }

17
libsrc/general/mysimd.hpp
View File

@ -54,8 +54,19 @@ namespace netgen
#ifdef __AVX__ #ifdef __AVX__
template <typename T>
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<> template<>
class alignas(32) SIMD<double> class alignas(32) SIMD<double> : public AlignedAlloc<SIMD<double>>
{ {
__m256d data; __m256d data;
@ -158,6 +169,10 @@ namespace netgen
#else #else
// it's only a dummy without AVX
template <typename T>
class AlignedAlloc { ; };
template<> template<>
class SIMD<double> class SIMD<double>
{ {

58
libsrc/general/ngpython.hpp
View File

@ -1,36 +1,42 @@
#ifdef NG_PYTHON #ifdef NG_PYTHON
#include <boost/python.hpp> #include <pybind11/pybind11.h>
namespace bp = boost::python; #include <pybind11/operators.h>
namespace py = pybind11;
#include <iostream> #include <iostream>
#include <sstream>
namespace pybind11 {
template<typename T>
bool CheckCast( py::handle obj ) {
try{
obj.cast<T>();
return true;
}
catch (py::cast_error &e) {
return false;
}
}
template <typename T>
struct extract
{
py::handle obj;
extract( py::handle aobj ) : obj(aobj) {}
bool check() { return CheckCast<T>(obj); }
T operator()() { return obj.cast<T>(); }
};
}
struct NGDummyArgument {};
inline void NOOP_Deleter(void *) { ; }
namespace netgen 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<std::string>(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 // Lambda to function pointer conversion
template <typename Function> template <typename Function>

7
libsrc/geom2d/geometry2d.hpp
View File

@ -109,7 +109,12 @@ namespace netgen
double curv = fabs(first(0)*second(1)-first(1)*second(0)) / pow(first.Length(), 3); double curv = fabs(first(0)*second(1)-first(1)*second(0)) / pow(first.Length(), 3);
return curv; return curv;
} }
virtual bool InConvexHull (Point<2> p, double eps) const
{
return seg.InConvexHull (p, eps);
}
}; };

129
libsrc/geom2d/python_geom2d.cpp
View File

@ -1,13 +1,11 @@
#ifdef NG_PYTHON #ifdef NG_PYTHON
#include <boost/python.hpp>
#include <../general/ngpython.hpp> #include <../general/ngpython.hpp>
#include <meshing.hpp> #include <meshing.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>
using namespace netgen; using namespace netgen;
namespace bp = boost::python;
namespace netgen namespace netgen
{ {
@ -15,24 +13,22 @@ namespace netgen
} }
DLL_HEADER void ExportGeom2d() DLL_HEADER void ExportGeom2d(py::module &m)
{ {
ModuleScope module("geom2d");
bp::class_<SplineGeometry2d, shared_ptr<SplineGeometry2d>, boost::noncopyable> py::class_<SplineGeometry2d, shared_ptr<SplineGeometry2d>>
("SplineGeometry", (m, "SplineGeometry",
"a 2d boundary representation geometry model by lines and splines") "a 2d boundary representation geometry model by lines and splines")
.def("__init__", bp::make_constructor .def(py::init<>())
(FunctionPointer .def("__init__",
([](const string & filename) [](SplineGeometry2d *instance, const string & filename)
{ {
cout << "load geometry"; cout << "load geometry";
ifstream ist(filename); ifstream ist(filename);
auto geom = make_shared<SplineGeometry2d>(); new (instance) SplineGeometry2d();
geom->Load (filename.c_str()); instance->Load (filename.c_str());
ng_geometry = geom; ng_geometry = shared_ptr<SplineGeometry2d>(instance, NOOP_Deleter);
return geom; })
})))
.def("Load",&SplineGeometry2d::Load) .def("Load",&SplineGeometry2d::Load)
.def("AppendPoint", FunctionPointer .def("AppendPoint", FunctionPointer
@ -47,16 +43,17 @@ DLL_HEADER void ExportGeom2d()
self.geompoints.Append(gp); self.geompoints.Append(gp);
return self.geompoints.Size()-1; return self.geompoints.Size()-1;
}), }),
(bp::arg("self"), bp::arg("x"), bp::arg("y"), bp::arg("maxh") = 1e99, bp::arg("hpref")=false)) py::arg("x"), py::arg("y"), py::arg("maxh") = 1e99, py::arg("hpref")=false)
.def("Append", FunctionPointer([](SplineGeometry2d &self, bp::list segment, int leftdomain, int rightdomain, bp::object bc, double maxh, bool hpref) .def("Append", FunctionPointer([](SplineGeometry2d &self, py::list segment, int leftdomain, int rightdomain,
py::object bc, py::object copy, double maxh, bool hpref)
{ {
bp::extract<std::string> segtype(segment[0]); py::extract<std::string> segtype(segment[0]);
SplineSegExt * seg; SplineSegExt * seg;
if (segtype().compare("line") == 0) if (segtype().compare("line") == 0)
{ {
bp::extract<int> point_index1(segment[1]); py::extract<int> point_index1(segment[1]);
bp::extract<int> point_index2(segment[2]); py::extract<int> point_index2(segment[2]);
//point_index1.check() //point_index1.check()
LineSeg<2> * l = new LineSeg<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2())); 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) else if (segtype().compare("spline3") == 0)
{ {
bp::extract<int> point_index1(segment[1]); py::extract<int> point_index1(segment[1]);
bp::extract<int> point_index2(segment[2]); py::extract<int> point_index2(segment[2]);
bp::extract<int> point_index3(segment[3]); py::extract<int> point_index3(segment[3]);
SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2()), self.GetPoint(point_index3())); SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1()), self.GetPoint(point_index2()), self.GetPoint(point_index3()));
seg = new SplineSegExt(*seg3); seg = new SplineSegExt(*seg3);
@ -82,11 +79,14 @@ DLL_HEADER void ExportGeom2d()
seg->hpref_right = hpref; seg->hpref_right = hpref;
seg->reffak = 1; seg->reffak = 1;
seg->copyfrom = -1; seg->copyfrom = -1;
if (bp::extract<int>(bc).check()) if (py::extract<int>(copy).check())
seg->bc = bp::extract<int>(bc)(); seg->copyfrom = py::extract<int>(copy)()+1;
else if (bp::extract<string>(bc).check())
if (py::extract<int>(bc).check())
seg->bc = py::extract<int>(bc)();
else if (py::extract<string>(bc).check())
{ {
string bcname = bp::extract<string>(bc)(); string bcname = py::extract<string>(bc)();
int bcnum = self.GetBCNumber(bcname); int bcnum = self.GetBCNumber(bcname);
if (bcnum == 0) if (bcnum == 0)
bcnum = self.AddBCName(bcname); bcnum = self.AddBCName(bcname);
@ -95,25 +95,26 @@ DLL_HEADER void ExportGeom2d()
else else
seg->bc = self.GetNSplines()+1; seg->bc = self.GetNSplines()+1;
self.AppendSegment(seg); self.AppendSegment(seg);
}), (bp::arg("self"), bp::arg("point_indices"), bp::arg("leftdomain") = 1, bp::arg("rightdomain") = 0, return self.GetNSplines()-1;
bp::arg("bc")=bp::object(), bp::arg("maxh")=1e99, bp::arg("hpref")=false }), 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; SplineSegExt * seg;
//int a = bp::extract<int>(point_indices[0]); //int a = py::extract<int>(point_indices[0]);
if (npts == 2) if (npts == 2)
{ {
LineSeg<2> * l = new LineSeg<2>(self.GetPoint(bp::extract<int>(point_indices[0])), self.GetPoint(bp::extract<int>(point_indices[1]))); LineSeg<2> * l = new LineSeg<2>(self.GetPoint(py::extract<int>(point_indices[0])()), self.GetPoint(py::extract<int>(point_indices[1])()));
seg = new SplineSegExt(*l); seg = new SplineSegExt(*l);
} }
else if (npts == 3) else if (npts == 3)
{ {
SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(bp::extract<int>(point_indices[0])), self.GetPoint(bp::extract<int>(point_indices[1])), self.GetPoint(bp::extract<int>(point_indices[2]))); SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(py::extract<int>(point_indices[0])()), self.GetPoint(py::extract<int>(point_indices[1])()), self.GetPoint(py::extract<int>(point_indices[2])()));
seg = new SplineSegExt(*seg3); seg = new SplineSegExt(*seg3);
} }
@ -123,7 +124,7 @@ DLL_HEADER void ExportGeom2d()
seg->reffak = 1; seg->reffak = 1;
seg->copyfrom = -1; seg->copyfrom = -1;
self.AppendSegment(seg); 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) //.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)); // 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; // seg->copyfrom = -1;
// self.AppendSegment(seg); // 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) //.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)); // 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->reffak = 1;
// seg->copyfrom = -1; // seg->copyfrom = -1;
// self.AppendSegment(seg); // 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) .def("SetMaterial", &SplineGeometry2d::SetMaterial)
@ -161,22 +162,22 @@ DLL_HEADER void ExportGeom2d()
Box<2> box(self.GetBoundingBox()); Box<2> box(self.GetBoundingBox());
double xdist = box.PMax()(0) - box.PMin()(0); double xdist = box.PMax()(0) - box.PMin()(0);
double ydist = box.PMax()(1) - box.PMin()(1); 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); py::tuple xlim = py::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 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++) 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) if (self.splines[i]->GetType().compare("line")==0)
{ {
GeomPoint<2> p1 = self.splines[i]->StartPI(); GeomPoint<2> p1 = self.splines[i]->StartPI();
GeomPoint<2> p2 = self.splines[i]->EndPI(); GeomPoint<2> p2 = self.splines[i]->EndPI();
xp.append(p1(0)); xp.append(py::cast(p1(0)));
xp.append(p2(0)); xp.append(py::cast(p2(0)));
yp.append(p1(1)); yp.append(py::cast(p1(1)));
yp.append(p2(1)); yp.append(py::cast(p2(1)));
} }
else if (self.splines[i]->GetType().compare("spline3")==0) else if (self.splines[i]->GetType().compare("spline3")==0)
{ {
@ -186,37 +187,37 @@ DLL_HEADER void ExportGeom2d()
for (int j = 0; j <= n; j++) for (int j = 0; j <= n; j++)
{ {
GeomPoint<2> point = self.splines[i]->GetPoint(j*1./n); GeomPoint<2> point = self.splines[i]->GetPoint(j*1./n);
xp.append(point(0)); xp.append(py::cast(point(0)));
yp.append(point(1)); yp.append(py::cast(point(1)));
} }
} }
else else
{ {
cout << "spline is neither line nor spline3" << endl; cout << "spline is neither line nor spline3" << endl;
} }
xpoints.append(xp); xpoints.append(py::cast(xp));
ypoints.append(yp); 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) .def("PointData", FunctionPointer([](SplineGeometry2d &self)
{ {
bp::list xpoints, ypoints, pointindex; py::list xpoints, ypoints, pointindex;
for (int i = 0; i < self.geompoints.Size(); i++) for (int i = 0; i < self.geompoints.Size(); i++)
{ {
pointindex.append(i); pointindex.append(py::cast(i));
xpoints.append(self.geompoints[i][0]); xpoints.append(py::cast(self.geompoints[i][0]));
ypoints.append(self.geompoints[i][1]); 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) .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++) for (int i = 0; i < self.splines.Size(); i++)
{ {
@ -226,13 +227,13 @@ DLL_HEADER void ExportGeom2d()
normal(0) = normal(1); normal(0) = normal(1);
normal(1) = -temp; normal(1) = -temp;
leftdom.append(self.GetSpline(i).leftdom); leftdom.append(py::cast(self.GetSpline(i).leftdom));
rightdom.append(self.GetSpline(i).rightdom); rightdom.append(py::cast(self.GetSpline(i).rightdom));
rightpoints.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(bp::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) .def("Print", FunctionPointer([](SplineGeometry2d &self)
@ -265,8 +266,10 @@ DLL_HEADER void ExportGeom2d()
} }
BOOST_PYTHON_MODULE(libgeom2d) { PYBIND11_PLUGIN(libgeom2d) {
ExportGeom2d(); py::module m("geom2d", "pybind geom2d");
ExportGeom2d(m);
return m.ptr();
} }
#endif #endif

5
libsrc/gprim/geomfuncs.hpp
View File

@ -82,9 +82,10 @@ namespace netgen
*/ */
// inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2) // inline Vec<3> Cross (const Vec<3> & v1, const Vec<3> & v2)
inline Vec<3> Cross (Vec<3> v1, Vec<3> v2) template <typename T>
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(1) * v2(2) - v1(2) * v2(1),
v1(2) * v2(0) - v1(0) * v2(2), v1(2) * v2(0) - v1(0) * v2(2),
v1(0) * v2(1) - v1(1) * v2(0) ); v1(0) * v2(1) - v1(1) * v2(0) );

34
libsrc/gprim/geomobjects.hpp
View File

@ -17,7 +17,7 @@ namespace netgen
template <int D, typename T> template <int D, typename T>
class Point class Point : public AlignedAlloc<Point<D,T>>
{ {
protected: protected:
@ -39,16 +39,18 @@ namespace netgen
Point (T ax, T ay, T az, T au) Point (T ax, T ay, T az, T au)
{ x[0] = ax; x[1] = ay; x[2] = az; x[3] = au;} { x[0] = ax; x[1] = ay; x[2] = az; x[3] = au;}
Point (const Point<D> & p2) template <typename T2>
{ for (int i = 0; i < D; i++) x[i] = p2.x[i]; } Point (const Point<D,T2> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2(i); }
explicit Point (const Vec<D> & v) explicit Point (const Vec<D> & v)
{ for (int i = 0; i < D; i++) x[i] = v(i); } { for (int i = 0; i < D; i++) x[i] = v(i); }
Point & operator= (const Point<D> & p2) template <typename T2>
Point & operator= (const Point<D,T2> & 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; return *this;
} }
@ -65,7 +67,7 @@ namespace netgen
}; };
template <int D, typename T> template <int D, typename T>
class Vec class Vec : public AlignedAlloc<Vec<D,T>>
{ {
protected: protected:
@ -90,17 +92,17 @@ namespace netgen
Vec (const Vec<D> & p2) Vec (const Vec<D> & p2)
{ for (int i = 0; i < D; i++) x[i] = p2.x[i]; } { for (int i = 0; i < D; i++) x[i] = p2.x[i]; }
explicit Vec (const Point<D> & p) explicit Vec (const Point<D,T> & p)
{ for (int i = 0; i < D; i++) x[i] = p(i); } { for (int i = 0; i < D; i++) x[i] = p(i); }
Vec (const Vec<D> & p1, const Vec<D> & p2) Vec (const Vec & p1, const Vec & p2)
{ for(int i=0; i<D; i++) x[i] = p2(i)-p1(1); } { for(int i=0; i<D; i++) x[i] = p2(i)-p1(1); }
template <typename T2>
Vec & operator= (const Vec<D> & p2) Vec & operator= (const Vec<D,T2> & 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; return *this;
} }
@ -131,12 +133,12 @@ namespace netgen
return l; return l;
} }
const Vec<D> & Normalize () Vec & Normalize ()
{ {
T l = Length(); T l = Length();
if (l != 0) // if (l != 0)
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
x[i] /= l; x[i] /= (l+1e-40);
return *this; return *this;
} }
@ -148,7 +150,7 @@ namespace netgen
template <int H, int W=H, typename T = double> template <int H, int W=H, typename T = double>
class Mat class Mat : public AlignedAlloc<Mat<H,W,T>>
{ {
protected: protected:

18
libsrc/gprim/geomops.hpp
View File

@ -29,10 +29,10 @@ namespace netgen
template <int D> template <int D, typename T>
inline Point<D> operator+ (const Point<D> & a, const Vec<D> & b) inline Point<D,T> operator+ (const Point<D,T> & a, const Vec<D,T> & b)
{ {
Point<D> res; Point<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) + b(i); res(i) = a(i) + b(i);
return res; return res;
@ -40,10 +40,10 @@ namespace netgen
template <int D> template <int D, typename T>
inline Vec<D> operator- (const Point<D> & a, const Point<D> & b) inline Vec<D,T> operator- (const Point<D,T> & a, const Point<D,T> & b)
{ {
Vec<D> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = a(i) - b(i); res(i) = a(i) - b(i);
return res; return res;
@ -69,10 +69,10 @@ namespace netgen
template <int D> template <int D, typename T>
inline Vec<D> operator* (double s, const Vec<D> & b) inline Vec<D,T> operator* (T s, const Vec<D,T> & b)
{ {
Vec<D> res; Vec<D,T> res;
for (int i = 0; i < D; i++) for (int i = 0; i < D; i++)
res(i) = s * b(i); res(i) = s * b(i);
return res; return res;

9
libsrc/gprim/geomtest3d.hpp
View File

@ -78,6 +78,15 @@ extern double MinDistLP2 (const Point3d & lp1, const Point3d & lp2, const Point3
extern double MinDistTP2 (const Point3d & tp1, const Point3d & tp2, extern double MinDistTP2 (const Point3d & tp1, const Point3d & tp2,
const Point3d & tp3, const Point3d & p); 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] /// Minimal distance of the 2 lines [l1p1,l1p2] and [l2p1,l2p2]
extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2, extern double MinDistLL2 (const Point3d & l1p1, const Point3d & l1p2,
const Point3d & l2p1, const Point3d & l2p2); const Point3d & l2p1, const Point3d & l2p2);

61
libsrc/gprim/spline.cpp
View File

@ -189,6 +189,51 @@ namespace netgen
if (tang * gradn < 0) u *= -1; if (tang * gradn < 0) u *= -1;
} }
template<int D>
void SplineSeg3<D> :: GetCoeff (Vector & u, Point<D> 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<D> 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<D> p0 = GetPoint(0);
Vec<D> 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; double t_old = -1;
/*
if(proj_latest_t > 0. && proj_latest_t < 1.) if(proj_latest_t > 0. && proj_latest_t < 1.)
t = proj_latest_t; t = proj_latest_t;
else else
t = 0.5; 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<D> phi; Point<D> phi;
Vec<D> phip,phipp,phimp; Vec<D> phip,phipp,phimp;

48
libsrc/gprim/spline.hpp
View File

@ -85,6 +85,7 @@ namespace netgen
void PrintCoeff (ostream & ost) const; void PrintCoeff (ostream & ost) const;
virtual void GetCoeff (Vector & coeffs) const = 0; virtual void GetCoeff (Vector & coeffs) const = 0;
virtual void GetCoeff (Vector & coeffs, Point<D> p0) const { ; }
virtual void GetPoints (int n, Array<Point<D> > & points) const; virtual void GetPoints (int n, Array<Point<D> > & points) const;
@ -95,6 +96,9 @@ namespace netgen
Array < Point<D> > & points, const double eps) const Array < Point<D> > & points, const double eps) const
{points.SetSize(0);} {points.SetSize(0);}
// is the point in the convex hull (increased by eps) of the spline ?
virtual bool InConvexHull (Point<D> p, double eps) const = 0;
virtual double MaxCurvature(void) const = 0; virtual double MaxCurvature(void) const = 0;
virtual string GetType(void) const {return "splinebase";} virtual string GetType(void) const {return "splinebase";}
@ -135,12 +139,18 @@ namespace netgen
virtual const GeomPoint<D> & EndPI () const { return p2; } virtual const GeomPoint<D> & EndPI () const { return p2; }
/// ///
virtual void GetCoeff (Vector & coeffs) const; virtual void GetCoeff (Vector & coeffs) const;
virtual void GetCoeff (Vector & coeffs, Point<D> p0) const;
virtual string GetType(void) const {return "line";} virtual string GetType(void) const {return "line";}
virtual void LineIntersections (const double a, const double b, const double c, virtual void LineIntersections (const double a, const double b, const double c,
Array < Point<D> > & points, const double eps) const; Array < Point<D> > & points, const double eps) const;
virtual bool InConvexHull (Point<D> p, double eps) const
{
return MinDistLP2 (p1, p2, p) < sqr(eps);
}
virtual double MaxCurvature(void) const {return 0;} virtual double MaxCurvature(void) const {return 0;}
virtual void Project (const Point<D> point, Point<D> & point_on_curve, double & t) const; virtual void Project (const Point<D> point, Point<D> & point_on_curve, double & t) const;
@ -178,7 +188,8 @@ namespace netgen
virtual const GeomPoint<D> & EndPI () const { return p3; } virtual const GeomPoint<D> & EndPI () const { return p3; }
/// ///
virtual void GetCoeff (Vector & coeffs) const; virtual void GetCoeff (Vector & coeffs) const;
virtual void GetCoeff (Vector & coeffs, Point<D> p0) const;
virtual string GetType(void) const {return "spline3";} virtual string GetType(void) const {return "spline3";}
const GeomPoint<D> & TangentPoint (void) const { return p2; } const GeomPoint<D> & 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, DLL_HEADER virtual void LineIntersections (const double a, const double b, const double c,
Array < Point<D> > & points, const double eps) const; Array < Point<D> > & points, const double eps) const;
virtual bool InConvexHull (Point<D> p, double eps) const
{
return MinDistTP2 (p1, p2, p3, p) < sqr(eps);
}
DLL_HEADER virtual double MaxCurvature(void) const; DLL_HEADER virtual double MaxCurvature(void) const;
DLL_HEADER virtual void Project (const Point<D> point, Point<D> & point_on_curve, double & t) const; DLL_HEADER virtual void Project (const Point<D> point, Point<D> & point_on_curve, double & t) const;
@ -232,6 +248,11 @@ namespace netgen
virtual void LineIntersections (const double a, const double b, const double c, virtual void LineIntersections (const double a, const double b, const double c,
Array < Point<D> > & points, const double eps) const; Array < Point<D> > & points, const double eps) const;
virtual bool InConvexHull (Point<D> p, double eps) const
{
return (Dist2 (p, pm) < sqr(radius+eps));
}
virtual double MaxCurvature(void) const {return 1./radius;} virtual double MaxCurvature(void) const {return 1./radius;}
}; };
@ -261,6 +282,10 @@ namespace netgen
virtual void GetCoeff (Vector & coeffs) const {;} virtual void GetCoeff (Vector & coeffs) const {;}
virtual double MaxCurvature(void) const {return 1;} virtual double MaxCurvature(void) const {return 1;}
// needs implementation ...
virtual bool InConvexHull (Point<D> p, double eps) const
{ return true; }
}; };
@ -372,6 +397,19 @@ namespace netgen
coeffs[5] = -dx * p1(1) + dy * p1(0); coeffs[5] = -dx * p1(1) + dy * p1(0);
} }
template<int D>
void LineSeg<D> :: GetCoeff (Vector & coeffs, Point<D> 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<int D> template<int D>
@ -597,6 +635,10 @@ namespace netgen
virtual void GetCoeff (Vector & coeffs) const {;} virtual void GetCoeff (Vector & coeffs) const {;}
virtual double MaxCurvature(void) const {return 1;} virtual double MaxCurvature(void) const {return 1;}
// needs implementation ...
virtual bool InConvexHull (Point<D> p, double eps) const
{ return true; }
}; };
// Constructor // Constructor

3
libsrc/gprim/transform3d.hpp
View File

@ -182,6 +182,9 @@ public:
{ {
to = m * from; to = m * from;
} }
Point<D> operator() (Point<D> from) const { Point<D> to; Transform(from, to); return to; }
Vec<D> operator() (Vec<D> from) const { Vec<D> to; Transform(from, to); return to; }
}; };
template <int D> template <int D>

4
libsrc/include/nginterface_v2_impl.hpp
View File

@ -64,6 +64,10 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<1> (int nr) const
Ng_Element ret; Ng_Element ret;
ret.type = NG_ELEMENT_TYPE(el.GetType()); ret.type = NG_ELEMENT_TYPE(el.GetType());
ret.index = el.si; 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.num = el.GetNP();
ret.points.ptr = (int*)&(el[0]); ret.points.ptr = (int*)&(el[0]);

2
libsrc/meshing/curvedelems.cpp
View File

@ -4176,6 +4176,7 @@ namespace netgen
T * dxdxi, size_t sdxdxi) T * dxdxi, size_t sdxdxi)
{ {
// multipointtrafovar++; // multipointtrafovar++;
/*
static int timer = NgProfiler::CreateTimer ("calcmultipointelementtrafo"); static int timer = NgProfiler::CreateTimer ("calcmultipointelementtrafo");
static int timer1 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 1"); static int timer1 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 1");
static int timer2 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 2"); static int timer2 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 2");
@ -4183,6 +4184,7 @@ namespace netgen
static int timer4 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 4"); static int timer4 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 4");
static int timer5 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 5"); static int timer5 = NgProfiler::CreateTimer ("calcmultipointelementtrafo 5");
NgProfiler::RegionTimer reg(timer); NgProfiler::RegionTimer reg(timer);
*/
// NgProfiler::StartTimer (timer); // NgProfiler::StartTimer (timer);
// NgProfiler::StartTimer (timer1); // NgProfiler::StartTimer (timer1);
if (mesh.coarsemesh) if (mesh.coarsemesh)

25
libsrc/meshing/meshtype.hpp
View File

@ -1020,13 +1020,6 @@ namespace netgen
class DLL_HEADER MeshingParameters class DLL_HEADER MeshingParameters
{ {
public: public:
/** /**
3d optimization strategy: 3d optimization strategy:
// m .. move nodes // m .. move nodes
@ -1130,12 +1123,28 @@ namespace netgen
MeshingParameters (); MeshingParameters ();
/// ///
MeshingParameters (const MeshingParameters & mp2) = default; 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 Print (ostream & ost) const;
/// ///
// void CopyFrom(const MeshingParameters & other); // 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<MeshSizePoint> meshsize_points;
void (*render_function)(bool) = NULL; void (*render_function)(bool) = NULL;
void Render(bool blocking = false) void Render(bool blocking = false)
{ {

433
libsrc/meshing/python_mesh.cpp
View File

@ -1,7 +1,5 @@
#ifdef NG_PYTHON #ifdef NG_PYTHON
#include <boost/python.hpp>
#include <boost/python/slice.hpp>
#include <../general/ngpython.hpp> #include <../general/ngpython.hpp>
#include <mystdlib.h> #include <mystdlib.h>
@ -9,10 +7,11 @@
#include <csg.hpp> #include <csg.hpp>
#include <geometry2d.hpp> #include <geometry2d.hpp>
#include <../interface/writeuser.hpp> #include <../interface/writeuser.hpp>
#include <nginterface.h>
#include <nginterface_v2.hpp>
using namespace netgen; using namespace netgen;
namespace bp = boost::python;
namespace netgen namespace netgen
{ {
@ -21,23 +20,19 @@ namespace netgen
template <typename T, int BASE = 0, typename TIND = int> template <typename T, int BASE = 0, typename TIND = int>
void ExportArray () void ExportArray (py::module &m)
{ {
string name = string("Array_") + typeid(T).name(); string name = string("Array_") + typeid(T).name();
bp::class_<Array<T,BASE,TIND>,boost::noncopyable>(name.c_str()) py::class_<Array<T,BASE,TIND>>(m, name.c_str())
.def ("__len__", &Array<T,BASE,TIND>::Size) .def ("__len__", [] ( Array<T,BASE,TIND> &self ) { return self.Size(); } )
.def ("__getitem__", .def ("__getitem__",
FunctionPointer ([](Array<T,BASE,TIND> & self, TIND i) -> T& FunctionPointer ([](Array<T,BASE,TIND> & self, TIND i) -> T&
{ {
if (i < BASE || i >= BASE+self.Size()) if (i < BASE || i >= BASE+self.Size())
bp::exec("raise IndexError()\n"); throw py::index_error();
return self[i]; return self[i];
}), }),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def ("__iter__",
bp::range (FunctionPointer([](Array<T,BASE,TIND> & self) { return &self[BASE]; }),
FunctionPointer([](Array<T,BASE,TIND> & self) { return &self[BASE+self.Size()]; })))
; ;
} }
@ -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<NgException>(&TranslateException);
bp::class_<PointIndex>("PointId", bp::init<int>()) py::class_<PointIndex>(m, "PointId")
.def(py::init<int>())
.def("__repr__", &ToString<PointIndex>) .def("__repr__", &ToString<PointIndex>)
.def("__str__", &ToString<PointIndex>) .def("__str__", &ToString<PointIndex>)
.add_property("nr", &PointIndex::operator int) .def_property_readonly("nr", &PointIndex::operator int)
.def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other) .def("__eq__" , FunctionPointer( [](PointIndex &self, PointIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) ) { return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast<int>(self); }) ) .def("__hash__" , FunctionPointer( [](PointIndex &self ) { return static_cast<int>(self); }) )
; ;
bp::class_<ElementIndex>("ElementId3D", bp::init<int>()) py::class_<ElementIndex>(m, "ElementId3D")
.def(py::init<int>())
.def("__repr__", &ToString<ElementIndex>) .def("__repr__", &ToString<ElementIndex>)
.def("__str__", &ToString<ElementIndex>) .def("__str__", &ToString<ElementIndex>)
.add_property("nr", &ElementIndex::operator int) .def_property_readonly("nr", &ElementIndex::operator int)
.def("__eq__" , FunctionPointer( [](ElementIndex &self, ElementIndex &other) .def("__eq__" , FunctionPointer( [](ElementIndex &self, ElementIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) ) { return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](ElementIndex &self ) { return static_cast<int>(self); }) ) .def("__hash__" , FunctionPointer( [](ElementIndex &self ) { return static_cast<int>(self); }) )
; ;
bp::class_<SurfaceElementIndex>("ElementId2D", bp::init<int>()) py::class_<SurfaceElementIndex>(m, "ElementId2D")
.def(py::init<int>())
.def("__repr__", &ToString<SurfaceElementIndex>) .def("__repr__", &ToString<SurfaceElementIndex>)
.def("__str__", &ToString<SurfaceElementIndex>) .def("__str__", &ToString<SurfaceElementIndex>)
.add_property("nr", &SurfaceElementIndex::operator int) .def_property_readonly("nr", &SurfaceElementIndex::operator int)
.def("__eq__" , FunctionPointer( [](SurfaceElementIndex &self, SurfaceElementIndex &other) .def("__eq__" , FunctionPointer( [](SurfaceElementIndex &self, SurfaceElementIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) ) { return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](SurfaceElementIndex &self ) { return static_cast<int>(self); }) ) .def("__hash__" , FunctionPointer( [](SurfaceElementIndex &self ) { return static_cast<int>(self); }) )
; ;
bp::class_<SegmentIndex>("ElementId1D", bp::init<int>()) py::class_<SegmentIndex>(m, "ElementId1D")
.def(py::init<int>())
.def("__repr__", &ToString<SegmentIndex>) .def("__repr__", &ToString<SegmentIndex>)
.def("__str__", &ToString<SegmentIndex>) .def("__str__", &ToString<SegmentIndex>)
.add_property("nr", &SegmentIndex::operator int) .def_property_readonly("nr", &SegmentIndex::operator int)
.def("__eq__" , FunctionPointer( [](SegmentIndex &self, SegmentIndex &other) .def("__eq__" , FunctionPointer( [](SegmentIndex &self, SegmentIndex &other)
{ return static_cast<int>(self)==static_cast<int>(other); }) ) { return static_cast<int>(self)==static_cast<int>(other); }) )
.def("__hash__" , FunctionPointer( [](SegmentIndex &self ) { return static_cast<int>(self); }) ) .def("__hash__" , FunctionPointer( [](SegmentIndex &self ) { return static_cast<int>(self); }) )
@ -96,176 +92,162 @@ DLL_HEADER void ExportNetgenMeshing()
/* /*
bp::class_<Point<3>> ("Point") py::class_<Point<3>> ("Point")
.def(bp::init<double,double,double>()) .def(py::init<double,double,double>())
; ;
*/ */
bp::class_<MeshPoint /* ,bp::bases<Point<3>> */ >("MeshPoint") py::class_<MeshPoint /* ,py::bases<Point<3>> */ >(m, "MeshPoint")
.def(bp::init<Point<3>>()) .def(py::init<Point<3>>())
.def("__str__", &ToString<MeshPoint>) .def("__str__", &ToString<MeshPoint>)
.def("__repr__", &ToString<MeshPoint>) .def("__repr__", &ToString<MeshPoint>)
.add_property("p", FunctionPointer([](const MeshPoint & self) .def_property_readonly("p", FunctionPointer([](const MeshPoint & self)
{ {
bp::list l; py::list l;
l.append ( self[0] ); l.append ( py::cast(self[0]) );
l.append ( self[1] ); l.append ( py::cast(self[1]) );
l.append ( self[2] ); l.append ( py::cast(self[2]) );
return bp::tuple(l); return py::tuple(l);
})) }))
.def("__getitem__", FunctionPointer([](const MeshPoint & self, int index) { .def("__getitem__", FunctionPointer([](const MeshPoint & self, int index) {
if(index<0 || index>2) if(index<0 || index>2)
bp::exec("raise IndexError()\n"); throw py::index_error();
return self[index]; return self[index];
})) }))
; ;
bp::class_<Element>("Element3D") py::class_<Element>(m, "Element3D")
.def("__init__", bp::make_constructor .def("__init__", [](Element *instance, int index, py::list vertices)
(FunctionPointer ([](int index, bp::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++) for (int i = 0; i < 4; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]); (*instance)[i] = py::extract<PointIndex>(vertices[i])();
tmp->SetIndex(index); instance->SetIndex(index);
return tmp;
} }
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++) for (int i = 0; i < 6; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]); (*instance)[i] = py::extract<PointIndex>(vertices[i])();
tmp->SetIndex(index); instance->SetIndex(index);
return tmp;
} }
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++) for (int i = 0; i < 8; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]); (*instance)[i] = py::extract<PointIndex>(vertices[i])();
tmp->SetIndex(index); instance->SetIndex(index);
return tmp;
} }
throw NgException ("cannot create element"); throw NgException ("cannot create element");
}), },
bp::default_call_policies(), // need it to use arguments py::arg("index")=1,py::arg("vertices"),
(bp::arg("index")=1,bp::arg("vertices"))),
"create volume element" "create volume element"
) )
.def("__repr__", &ToString<Element>) .def("__repr__", &ToString<Element>)
.add_property("index", &Element::GetIndex, &Element::SetIndex) .def_property("index", &Element::GetIndex, &Element::SetIndex)
.add_property("vertices", .def_property_readonly("vertices",
FunctionPointer ([](const Element & self) -> bp::list FunctionPointer ([](const Element & self) -> py::list
{ {
bp::list li; py::list li;
for (int i = 0; i < self.GetNV(); i++) for (int i = 0; i < self.GetNV(); i++)
li.append (self[i]); li.append (py::cast(self[i]));
return li; return li;
})) }))
; ;
bp::class_<Element2d>("Element2D") py::class_<Element2d>(m, "Element2D")
.def("__init__", bp::make_constructor .def("__init__",
(FunctionPointer ([](int index, bp::list vertices) [](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++) for (int i = 0; i < 3; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]); (*instance)[i] = py::extract<PointIndex>(vertices[i])();
tmp->SetIndex(index); instance->SetIndex(index);
return tmp;
} }
else else
{ {
Element2d * tmp = new Element2d(QUAD); new (instance) Element2d(QUAD);
for (int i = 0; i < 4; i++) for (int i = 0; i < 4; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]); (*instance)[i] = py::extract<PointIndex>(vertices[i])();
tmp->SetIndex(index); instance->SetIndex(index);
return tmp;
} }
}), },
bp::default_call_policies(), // need it to use arguments py::arg("index")=1,py::arg("vertices"),
(bp::arg("index")=1,bp::arg("vertices"))),
"create surface element" "create surface element"
) )
.add_property("index", &Element2d::GetIndex, &Element2d::SetIndex) .def_property("index", &Element2d::GetIndex, &Element2d::SetIndex)
.add_property("vertices", .def_property_readonly("vertices",
FunctionPointer([](const Element2d & self) -> bp::list FunctionPointer([](const Element2d & self) -> py::list
{ {
bp::list li; py::list li;
for (int i = 0; i < self.GetNV(); i++) for (int i = 0; i < self.GetNV(); i++)
li.append(self[i]); li.append(py::cast(self[i]));
return li; return li;
})) }))
; ;
bp::class_<Segment>("Element1D") py::class_<Segment>(m, "Element1D")
.def("__init__", bp::make_constructor .def("__init__",
(FunctionPointer ([](bp::list vertices, bp::list surfaces, int index) [](Segment *instance, py::list vertices, py::list surfaces, int index)
{ {
Segment * tmp = new Segment; new (instance) Segment();
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
(*tmp)[i] = bp::extract<PointIndex>(vertices[i]); (*instance)[i] = py::extract<PointIndex>(vertices[i])();
tmp -> si = index; instance -> si = index;
if (len(surfaces)) if (len(surfaces))
{ {
tmp->surfnr1 = bp::extract<int>(surfaces[0]); instance->surfnr1 = py::extract<int>(surfaces[0])();
tmp->surfnr2 = bp::extract<int>(surfaces[1]); instance->surfnr2 = py::extract<int>(surfaces[1])();
} }
return tmp; },
}), py::arg("vertices"),
bp::default_call_policies(), py::arg("surfaces")=py::list(),
(bp::arg("vertices"), py::arg("index")=1,
bp::arg("surfaces")=bp::list(),
bp::arg("index")=1
)),
"create segment element" "create segment element"
) )
.def("__repr__", &ToString<Segment>) .def("__repr__", &ToString<Segment>)
.add_property("vertices", .def_property_readonly("vertices",
FunctionPointer ([](const Segment & self) -> bp::list FunctionPointer ([](const Segment & self) -> py::list
{ {
bp::list li; py::list li;
for (int i = 0; i < 2; i++) for (int i = 0; i < 2; i++)
li.append (self[i]); li.append (py::cast(self[i]));
return li; return li;
})) }))
.add_property("surfaces", .def_property_readonly("surfaces",
FunctionPointer ([](const Segment & self) -> bp::list FunctionPointer ([](const Segment & self) -> py::list
{ {
bp::list li; py::list li;
li.append (self.surfnr1); li.append (py::cast(self.surfnr1));
li.append (self.surfnr2); li.append (py::cast(self.surfnr2));
return li; return li;
})) }))
; ;
bp::class_<Element0d>("Element0D") py::class_<Element0d>(m, "Element0D")
.def("__init__", bp::make_constructor .def("__init__",
(FunctionPointer ([](PointIndex vertex, int index) [](Element0d *instance, PointIndex vertex, int index)
{ {
Element0d * tmp = new Element0d; new (instance) Element0d;
tmp->pnum = vertex; instance->pnum = vertex;
tmp->index = index; instance->index = index;
return tmp; },
}), py::arg("vertex"),
bp::default_call_policies(), py::arg("index")=1,
(bp::arg("vertex"),
bp::arg("index")=1
)),
"create point element" "create point element"
) )
.def("__repr__", &ToString<Element0d>) .def("__repr__", &ToString<Element0d>)
.add_property("vertices", .def_property_readonly("vertices",
FunctionPointer ([](const Element0d & self) -> bp::list FunctionPointer ([](const Element0d & self) -> py::list
{ {
bp::list li; py::list li;
li.append (self.pnum); li.append (py::cast(self.pnum));
return li; return li;
})) }))
; ;
@ -274,62 +256,65 @@ DLL_HEADER void ExportNetgenMeshing()
bp::class_<FaceDescriptor>("FaceDescriptor") py::class_<FaceDescriptor>(m, "FaceDescriptor")
.def(bp::init<const FaceDescriptor&>()) .def(py::init<const FaceDescriptor&>())
.def("__init__", bp::make_constructor .def("__init__",
(FunctionPointer ([](int surfnr, int domin, int domout, int bc) [](FaceDescriptor *instance, int surfnr, int domin, int domout, int bc)
{ {
auto fd = new FaceDescriptor(); new (instance) FaceDescriptor();
fd->SetSurfNr(surfnr); instance->SetSurfNr(surfnr);
fd->SetDomainIn(domin); instance->SetDomainIn(domin);
fd->SetDomainOut(domout); instance->SetDomainOut(domout);
fd->SetBCProperty(bc); instance->SetBCProperty(bc);
return fd; },
}), py::arg("surfnr")=1,
bp::default_call_policies(), // need it to use arguments py::arg("domin")=1,
(bp::arg("surfnr")=1, py::arg("domout")=py::int_(0),
bp::arg("domin")=1, py::arg("bc")=py::int_(0),
bp::arg("domout")=0,
bp::arg("bc")=0
)),
"create facedescriptor") "create facedescriptor")
.def("__str__", &ToString<FaceDescriptor>) .def("__str__", &ToString<FaceDescriptor>)
.def("__repr__", &ToString<FaceDescriptor>) .def("__repr__", &ToString<FaceDescriptor>)
.add_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr) .def_property("surfnr", &FaceDescriptor::SurfNr, &FaceDescriptor::SetSurfNr)
.add_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn) .def_property("domin", &FaceDescriptor::DomainIn, &FaceDescriptor::SetDomainIn)
.add_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut) .def_property("domout", &FaceDescriptor::DomainOut, &FaceDescriptor::SetDomainOut)
.add_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty) .def_property("bc", &FaceDescriptor::BCProperty, &FaceDescriptor::SetBCProperty)
.add_property("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); })) .def_property_readonly("bcname", FunctionPointer ([](FaceDescriptor & self) -> string { return self.GetBCName(); }))
.def("SetSurfaceColor", [](FaceDescriptor & self, py::list color )
{
Vec3d c;
c.X() = py::extract<double>(color[0])();
c.Y() = py::extract<double>(color[1])();
c.Z() = py::extract<double>(color[2])();
self.SetSurfColour(c);
})
; ;
ExportArray<Element>(); ExportArray<Element>(m);
ExportArray<Element2d>(); ExportArray<Element2d>(m);
ExportArray<Segment>(); ExportArray<Segment>(m);
ExportArray<Element0d>(); ExportArray<Element0d>(m);
ExportArray<MeshPoint,PointIndex::BASE,PointIndex>(); ExportArray<MeshPoint,PointIndex::BASE,PointIndex>(m);
ExportArray<FaceDescriptor>(); ExportArray<FaceDescriptor>(m);
;
#if (BOOST_VERSION >= 106000) && (BOOST_VERSION < 106100)
bp::register_ptr_to_python<shared_ptr<Mesh>>();
#endif
bp::class_<Mesh,shared_ptr<Mesh>,boost::noncopyable>("Mesh", bp::no_init)
// .def(bp::init<>("create empty mesh"))
.def("__init__", bp::make_constructor py::implicitly_convertible< int, PointIndex>();
(FunctionPointer ([](int dim)
py::class_<Ngx_Mesh>(m, "Ngx_Mesh")
.def(py::init<shared_ptr<Mesh>>())
.def_property_readonly("ngmesh", &Ngx_Mesh::GetMesh)
;
py::class_<Mesh,shared_ptr<Mesh>>(m, "Mesh")
// .def(py::init<>("create empty mesh"))
.def("__init__",
[](Mesh *instance, int dim)
{ {
auto mesh = make_shared<Mesh>(); new (instance) Mesh();
mesh->SetDimension(dim); instance->SetDimension(dim);
return mesh; },
}), py::arg("dim")=3
bp::default_call_policies(), // need it to use named arguments )
(
bp::arg("dim")=3
)
))
.def("__str__", &ToString<Mesh>) .def("__str__", &ToString<Mesh>)
@ -371,8 +356,8 @@ DLL_HEADER void ExportNetgenMeshing()
})) }))
// static_cast<void(Mesh::*)(const string & name)>(&Mesh::Load)) // static_cast<void(Mesh::*)(const string & name)>(&Mesh::Load))
.def("Save", static_cast<void(Mesh::*)(const string & name)const>(&Mesh::Save)) .def("Save", static_cast<void(Mesh::*)(const string & name)const>(&Mesh::Save))
.def("Export", FunctionPointer .def("Export",
([] (Mesh & self, string filename, string format) [] (Mesh & self, string filename, string format)
{ {
if (WriteUserFormat (format, self, *self.GetGeometry(), filename)) if (WriteUserFormat (format, self, *self.GetGeometry(), filename))
{ {
@ -384,35 +369,35 @@ DLL_HEADER void ExportNetgenMeshing()
err += string("'") + name + "'\n"; err += string("'") + name + "'\n";
throw NgException (err); 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", .def("Elements3D",
static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements), static_cast<Array<Element>&(Mesh::*)()> (&Mesh::VolumeElements),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def("Elements2D", .def("Elements2D",
static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements), static_cast<Array<Element2d>&(Mesh::*)()> (&Mesh::SurfaceElements),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def("Elements1D", .def("Elements1D",
static_cast<Array<Segment>&(Mesh::*)()> (&Mesh::LineSegments), static_cast<Array<Segment>&(Mesh::*)()> (&Mesh::LineSegments),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def("Elements0D", FunctionPointer([] (Mesh & self) -> Array<Element0d>& .def("Elements0D", FunctionPointer([] (Mesh & self) -> Array<Element0d>&
{ {
return self.pointelements; return self.pointelements;
} ), } ),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def("Points", .def("Points",
static_cast<Mesh::T_POINTS&(Mesh::*)()> (&Mesh::Points), static_cast<Mesh::T_POINTS&(Mesh::*)()> (&Mesh::Points),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor), .def("FaceDescriptor", static_cast<FaceDescriptor&(Mesh::*)(int)> (&Mesh::GetFaceDescriptor),
bp::return_value_policy<bp::reference_existing_object>()) py::return_value_policy::reference)
.def("GetNFaceDescriptors", &Mesh::GetNFD) .def("GetNFaceDescriptors", &Mesh::GetNFD)
@ -457,16 +442,23 @@ DLL_HEADER void ExportNetgenMeshing()
.def ("SetMaterial", &Mesh::SetMaterial) .def ("SetMaterial", &Mesh::SetMaterial)
.def ("GetMaterial", FunctionPointer([](Mesh & self, int domnr) .def ("GetMaterial", FunctionPointer([](Mesh & self, int domnr)
{ return string(self.GetMaterial(domnr)); })) { return string(self.GetMaterial(domnr)); }))
.def ("GenerateVolumeMesh", FunctionPointer .def ("GenerateVolumeMesh",
([](Mesh & self) [](Mesh & self, py::object pymp)
{ {
cout << "generate vol mesh" << endl; cout << "generate vol mesh" << endl;
MeshingParameters mp; MeshingParameters mp;
mp.optsteps3d = 5; if (py::extract<MeshingParameters>(pymp).check())
mp = py::extract<MeshingParameters>(pymp)();
else
{
mp.optsteps3d = 5;
}
MeshVolume (mp, self); MeshVolume (mp, self);
OptimizeVolume (mp, self); OptimizeVolume (mp, self);
})) },
py::arg("mp")=NGDummyArgument())
.def ("OptimizeVolumeMesh", FunctionPointer .def ("OptimizeVolumeMesh", FunctionPointer
([](Mesh & self) ([](Mesh & self)
@ -491,18 +483,19 @@ DLL_HEADER void ExportNetgenMeshing()
self.SetGeometry(geo); self.SetGeometry(geo);
})) }))
.def ("SetGeometry", FunctionPointer // TODO: fix this dependency on libgeom2d.so
([](Mesh & self, shared_ptr<SplineGeometry2d> geo) // .def ("SetGeometry", FunctionPointer
{ // ([](Mesh & self, shared_ptr<SplineGeometry2d> geo)
self.SetGeometry(geo); // {
})) // self.SetGeometry(geo);
// }))
.def ("BuildSearchTree", &Mesh::BuildElementSearchTree) .def ("BuildSearchTree", &Mesh::BuildElementSearchTree)
.def ("BoundaryLayer", FunctionPointer .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; BoundaryLayerParameters blp;
for (int i = 1; i <= self.GetNFD(); i++) for (int i = 1; i <= self.GetNFD(); i++)
@ -521,9 +514,9 @@ DLL_HEADER void ExportNetgenMeshing()
cout << "maxind = " << maxind << endl; cout << "maxind = " << maxind << endl;
for ( int i=0; i<n; i++ ) for ( int i=0; i<n; i++ )
{ {
blp.heights.Append( bp::extract<double>(thicknesses[i])()) ; blp.heights.Append( py::extract<double>(thicknesses[i])()) ;
blp.new_matnrs.Append( maxind+1+i ); blp.new_matnrs.Append( maxind+1+i );
self.SetMaterial (maxind+1+i, bp::extract<string>(materials[i])().c_str()); self.SetMaterial (maxind+1+i, py::extract<string>(materials[i])().c_str());
} }
blp.bulk_matnr = volnr; blp.bulk_matnr = volnr;
GenerateBoundaryLayer (self, blp); GenerateBoundaryLayer (self, blp);
@ -558,39 +551,48 @@ DLL_HEADER void ExportNetgenMeshing()
)) ))
; ;
py::implicitly_convertible< shared_ptr<Mesh>, Ngx_Mesh >();
typedef MeshingParameters MP; typedef MeshingParameters MP;
bp::class_<MP> ("MeshingParameters", bp::init<>()) py::class_<MP> (m, "MeshingParameters")
.def("__init__", bp::make_constructor .def(py::init<>())
(FunctionPointer ([](double maxh, bool quad_dominated) .def("__init__",
[](MP *instance, double maxh, bool quad_dominated, int optsteps2d, int optsteps3d)
{ {
auto tmp = new MeshingParameters; new (instance) MeshingParameters;
tmp->maxh = maxh; instance->maxh = maxh;
tmp->quad = int(quad_dominated); instance->quad = int(quad_dominated);
return tmp; instance->optsteps2d = optsteps2d;
}), instance->optsteps3d = optsteps3d;
bp::default_call_policies(), // need it to use arguments },
( py::arg("maxh")=1000,
bp::arg("maxh")=1000, py::arg("quad_dominated")=false,
bp::arg("quad_dominated")=false py::arg("optsteps2d") = 3,
)), py::arg("optsteps3d") = 3
,
"create meshing parameters" "create meshing parameters"
) )
.def("__str__", &ToString<MP>) .def("__str__", &ToString<MP>)
.add_property("maxh", .def_property("maxh",
FunctionPointer ([](const MP & mp ) { return mp.maxh; }), FunctionPointer ([](const MP & mp ) { return mp.maxh; }),
FunctionPointer ([](MP & mp, double maxh) { return mp.maxh = 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; delete testout;
testout = new ofstream (filename); testout = new ofstream (filename);
})); }));
bp::def("SetMessageImportance", FunctionPointer ([] (int importance) m.def("SetMessageImportance", FunctionPointer ([] (int importance)
{ {
int old = printmessage_importance; int old = printmessage_importance;
printmessage_importance = importance; printmessage_importance = importance;
@ -598,14 +600,11 @@ DLL_HEADER void ExportNetgenMeshing()
})); }));
} }
PYBIND11_PLUGIN(libmesh) {
py::module m("mesh", "pybind mesh");
BOOST_PYTHON_MODULE(libmesh) { ExportNetgenMeshing(m);
ExportNetgenMeshing(); return m.ptr();
} }
#endif #endif

357
libsrc/meshing/topology.cpp
View File

@ -62,7 +62,63 @@ namespace netgen
delete vert2pointelement; delete vert2pointelement;
} }
template <typename FUNC>
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 <typename FUNC> template <typename FUNC>
void LoopOverFaces (const Mesh & mesh, MeshTopology & top, PointIndex v, void LoopOverFaces (const Mesh & mesh, MeshTopology & top, PointIndex v,
FUNC func) FUNC func)
@ -139,11 +195,12 @@ namespace netgen
swap (face4.I2(), face4.I4()); swap (face4.I2(), face4.I4());
} }
if (face4.I1() != v) continue;
func(face4, elnr, j, true, facedir); func(face4, elnr, j, true, facedir);
/* /*
INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); INDEX_3 face(face4.I1(), face4.I2(), face4.I3());
if (face.I1() != v) continue;
if (vert2face.Used (face)) if (vert2face.Used (face))
{ {
@ -256,10 +313,10 @@ namespace netgen
swap (face4.I2(), face4.I4()); swap (face4.I2(), face4.I4());
} }
if (face4.I1() != v) continue;
func(face4, elnr, 0, false, facedir); func(face4, elnr, 0, false, facedir);
/* /*
INDEX_3 face(face4.I1(), face4.I2(), face4.I3()); INDEX_3 face(face4.I1(), face4.I2(), face4.I3());
if (face.I1() != v) continue;
if (vert2face.Used (face)) if (vert2face.Used (face))
facenum = vert2face.Get(face); facenum = vert2face.Get(face);
@ -300,7 +357,6 @@ namespace netgen
int nseg = mesh.GetNSeg(); int nseg = mesh.GetNSeg();
int np = mesh.GetNP(); int np = mesh.GetNP();
int nv = mesh.GetNV(); int nv = mesh.GetNV();
int ned = edge2vert.Size();
if (id == 0) if (id == 0)
PrintMessage (3, "Update mesh topology"); PrintMessage (3, "Update mesh topology");
@ -417,7 +473,7 @@ namespace netgen
cnt[edge2vert[i][0]]++; cnt[edge2vert[i][0]]++;
TABLE<int,PointIndex::BASE> vert2edge (cnt); TABLE<int,PointIndex::BASE> vert2edge (cnt);
for (int i = 0; i < edge2vert.Size(); i++) 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 // ensure all coarse grid and intermediate level edges
cnt = 0; cnt = 0;
@ -433,167 +489,141 @@ namespace netgen
if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]); if (parents[0] > PointIndex::BASE) vert2vertcoarse.AddSave (parents[0], parents[1]);
} }
Array<int,PointIndex::BASE> edgenr(nv);
Array<int,PointIndex::BASE> edgeflag(nv);
Array<int> vertex2;
edgeflag = PointIndex::BASE-1;
ned = edge2vert.Size();
int max_edge_on_vertex = 0;
for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++) for (int i = PointIndex::BASE; i < nv+PointIndex::BASE; i++)
{ {
vertex2.SetSize (0); int onv = vert2edge[i].Size() + vert2vertcoarse[i].Size() +
4*(*vert2element)[i].Size() + 2*(*vert2surfelement)[i].Size() + (*vert2segment)[i].Size();
for (int j = 0; j < vert2edge[i].Size(); j++) max_edge_on_vertex = max (onv, max_edge_on_vertex);
{
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<ElementIndex> 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;
}
} }
}
// count edges associated with vertices
cnt = 0;
ParallelForRange
(tm, mesh.Points().Size(),
[&] (size_t begin, size_t end)
{
INDEX_CLOSED_HASHTABLE<int> 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<int> v2eht(2*max_edge_on_vertex+10);
// Array<int> 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<int> v2eht(2*max_edge_on_vertex+10);
Array<int> 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 // generate faces
@ -638,7 +668,6 @@ namespace netgen
NgProfiler::StopTimer (timer2a); NgProfiler::StopTimer (timer2a);
NgProfiler::StartTimer (timer2b); NgProfiler::StartTimer (timer2b);
@ -764,8 +793,9 @@ namespace netgen
}); });
} }
}); });
/*
int oldnfa = face2vert.Size(); /*
int oldnfa = face2vert.Size();
int nfa = oldnfa; int nfa = oldnfa;
INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10); INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*max_face_on_vertex+10);
@ -1010,9 +1040,8 @@ namespace netgen
} }
} }
} }
*/
face2vert.SetAllocSize (nfa); face2vert.SetAllocSize (nfa);
*/
// *testout << "face2vert = " << endl << face2vert << endl; // *testout << "face2vert = " << endl << face2vert << endl;

5
libsrc/stlgeom/stltool.cpp
View File

@ -615,7 +615,7 @@ STLChart :: STLChart(STLGeometry * ageometry)
geometry = ageometry; geometry = ageometry;
if ( (stlparam.usesearchtree == 1)) if ( stlparam.usesearchtree == 1)
searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1), searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
geometry->GetBoundingBox().PMax() + Vec3d(1,1,1)); geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
else 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++) for (i = 0; i <= divisions; i++)
{ {

5
libsrc/visualization/mvdraw.cpp
View File

@ -776,13 +776,16 @@ namespace netgen
{ {
string displname; string displname;
// Display * dpy = glXGetCurrentDisplay(); Display * dpy = glXGetCurrentDisplay();
GLXDrawable drawable = glXGetCurrentDrawable(); GLXDrawable drawable = glXGetCurrentDrawable();
GLXContext ctx = glXGetCurrentContext(); GLXContext ctx = glXGetCurrentContext();
GLXContextID xid = glXGetContextIDEXT (ctx); GLXContextID xid = glXGetContextIDEXT (ctx);
displname = XDisplayName (0); 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 << "Init Parallel GL" << endl;
cout << "DisplayName = " << displname << endl; cout << "DisplayName = " << displname << endl;

12
libsrc/visualization/soldata.hpp
View File

@ -92,6 +92,18 @@ namespace netgen
return res; 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) virtual bool GetSegmentValue (int segnr, double xref, double * values)
{ return false; } { return false; }

29
libsrc/visualization/vsmesh.cpp
View File

@ -3529,14 +3529,11 @@ namespace netgen
#ifdef NG_PYTHON #ifdef NG_PYTHON
#include <boost/python.hpp>
#include <../general/ngpython.hpp> #include <../general/ngpython.hpp>
namespace bp = boost::python;
DLL_HEADER void ExportMeshVis() DLL_HEADER void ExportMeshVis(py::module &m)
{ {
using namespace netgen; using namespace netgen;
ModuleScope myscope("meshvis");
vispar.drawcolorbar = true; vispar.drawcolorbar = true;
vispar.drawnetgenlogo = true; vispar.drawnetgenlogo = true;
vispar.drawcoordinatecross = true; vispar.drawcoordinatecross = true;
@ -3546,12 +3543,12 @@ DLL_HEADER void ExportMeshVis()
vispar.drawtets = true; vispar.drawtets = true;
vispar.drawprisms = true; vispar.drawprisms = true;
vispar.drawoutline = true; vispar.drawoutline = true;
bp::class_<VisualSceneMesh, shared_ptr<VisualSceneMesh>> py::class_<VisualSceneMesh, shared_ptr<VisualSceneMesh>>
("VisualSceneMesh", bp::no_init) (m, "VisualSceneMesh")
.def("Draw", &VisualSceneMesh::DrawScene) .def("Draw", &VisualSceneMesh::DrawScene)
; ;
bp::def("VS", FunctionPointer m.def("VS", FunctionPointer
([](shared_ptr<Mesh> mesh) ([](shared_ptr<Mesh> mesh)
{ {
auto vs = make_shared<VisualSceneMesh>(); auto vs = make_shared<VisualSceneMesh>();
@ -3560,15 +3557,23 @@ DLL_HEADER void ExportMeshVis()
return vs; return vs;
})); }));
bp::def("MouseMove", FunctionPointer m.def("MouseMove", FunctionPointer
([](VisualSceneMesh &vsmesh, int oldx, int oldy, int newx, int ([](VisualSceneMesh &vsmesh, int oldx, int oldy, int newx, int
newy, char mode) newy, char mode)
{ {
vsmesh.MouseMove(oldx, oldy, newx, newy, 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) // BOOST_PYTHON_MODULE(libvisual)
{ // {
ExportMeshVis(); // ExportMeshVis();
} // }
#endif #endif

378
libsrc/visualization/vssolution.cpp
View File

@ -8,10 +8,7 @@
// #include <parallel.hpp> // #include <parallel.hpp>
#include <visual.hpp> #include <visual.hpp>
#include <limits> #include <limits>
namespace netgen namespace netgen
{ {
@ -441,6 +438,20 @@ namespace netgen
// glEnable(GL_BLEND); // glEnable(GL_BLEND);
glDisable(GL_BLEND); glDisable(GL_BLEND);
glCallList (surfellist); 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); glDisable(GL_BLEND);
/* /*
// transparent test ... // transparent test ...
@ -1182,6 +1193,17 @@ namespace netgen
shared_ptr<Mesh> mesh = GetMesh(); shared_ptr<Mesh> mesh = GetMesh();
static int timer = NgProfiler::CreateTimer ("Solution::DrawSurfaceElements"); 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); NgProfiler::RegionTimer reg (timer);
@ -1213,7 +1235,7 @@ namespace netgen
} }
#endif #endif
NgProfiler::StartTimer(timerstart);
if (surfellist) if (surfellist)
glDeleteLists (surfellist, 1); glDeleteLists (surfellist, 1);
@ -1256,10 +1278,58 @@ namespace netgen
Array<double> values(npt); Array<double> values(npt);
Array<double> mvalues(npt); Array<double> mvalues(npt);
int sol_comp = (sol && sol->draw_surface) ? sol->components : 0;
#ifdef __AVX__
Array<Point<2,SIMD<double>> > simd_pref ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<Point<3,SIMD<double>> > simd_points ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<Mat<3,2,SIMD<double>> > simd_dxdxis ( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<Vec<3,SIMD<double>> > simd_nvs( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() );
Array<SIMD<double>> simd_values( (npt+SIMD<double>::Size()-1)/SIMD<double>::Size() * sol_comp);
#endif
// Array<Point<3,float>> glob_pnts;
// Array<Vec<3,float>> glob_nvs;
// Array<double> glob_values;
if (sol && sol->draw_surface) mvalues.SetSize (npt * sol->components); if (sol && sol->draw_surface) mvalues.SetSize (npt * sol->components);
Array<complex<double> > valuesc(npt); Array<complex<double> > 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++) for (SurfaceElementIndex sei = 0; sei < nse; sei++)
{ {
const Element2d & el = (*mesh)[sei]; const Element2d & el = (*mesh)[sei];
@ -1413,7 +1483,38 @@ namespace netgen
n = 1 << subdivisions; n = 1 << subdivisions;
double invn = 1.0 / n; double invn = 1.0 / n;
npt = (n+1)*(n+2)/2; 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<double>::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<int> 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<int> glob_ind;
glob_ind.SetSize(ind_reftrig.Size());
for(SurfaceElementIndex sei = 0; sei < nse; sei++) for(SurfaceElementIndex sei = 0; sei < nse; sei++)
{ {
const Element2d & el = (*mesh)[sei]; const Element2d & el = (*mesh)[sei];
@ -1436,6 +1537,90 @@ namespace netgen
if ( el.GetType() == TRIG || el.GetType() == TRIG6 ) 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<double>> p1 = mesh->Point (el[0]);
Point<3,SIMD<double>> p2 = mesh->Point (el[1]);
Point<3,SIMD<double>> p3 = mesh->Point (el[2]);
Vec<3,SIMD<double>> vx = p1-p3;
Vec<3,SIMD<double>> 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<double>> 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); bool curved = curv.IsSurfaceElementCurved(sei);
for (int iy = 0, ii = 0; iy <= n; iy++) for (int iy = 0, ii = 0; iy <= n; iy++)
@ -1492,6 +1677,8 @@ namespace netgen
if (deform) if (deform)
for (int ii = 0; ii < npt; ii++) for (int ii = 0; ii < npt; ii++)
points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1)); points[ii] += GetSurfDeformation (sei, -1, pref[ii](0), pref[ii](1));
#endif
NgProfiler::StopTimer(timer1);
int save_usetexture = usetexture; int save_usetexture = usetexture;
if (!drawelem) if (!drawelem)
@ -1500,108 +1687,62 @@ namespace netgen
SetTextureMode (usetexture); 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); glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[0]);
for (int ix = 0; ix <= n-iy; ix++, ii++) glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(Point<3,double>),
for (int k = 0; k < 2; k++) npt*sizeof(Point<3,double>), &points[0][0]);
{ glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[1]);
if (ix+iy+k > n) continue; glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(Vec<3,double>),
int hi = (k == 0) ? ii : ii+n-iy+1; npt*sizeof(Vec<3,double>), &nvs[0][0]);
glBindBuffer (GL_ARRAY_BUFFER, surfel_vbo[2]);
if (drawelem) glBufferSubData (GL_ARRAY_BUFFER, base_pi*sizeof(double),
{ npt*sizeof(double), &values[0]);
if (usetexture != 2)
SetOpenGlColor (values[hi]); for (size_t i = 0; i < ind_reftrig.Size(); i++)
else glob_ind[i] = base_pi+ind_reftrig[i];
glTexCoord2f ( valuesc[hi].real(), valuesc[hi].imag() );
} glBindBuffer (GL_ELEMENT_ARRAY_BUFFER, surfel_vbo[3]);
else glBufferSubData (GL_ELEMENT_ARRAY_BUFFER, surfel_vbo_size*sizeof(int),
glColor4fv (col_grey); ind_reftrig.Size()*sizeof(int), &glob_ind[0]);
surfel_vbo_size += ind_reftrig.Size();
glNormal3dv (nvs[hi]); base_pi += npt;
glVertex3dv (points[hi]);
}
glEnd();
} }
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<int> 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)) if (!drawelem && (usetexture != save_usetexture))
{ {
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 (); glEndList ();
NgProfiler::StopTimer(timerlist);
#ifdef PARALLELGL #ifdef PARALLELGL
glFinish(); glFinish();
if (id > 0) if (id > 0)
@ -2492,12 +2656,14 @@ namespace netgen
if(minv_local < minv) if(minv_local < minv)
{ {
lock_guard<mutex> guard(min_mutex); lock_guard<mutex> guard(min_mutex);
minv = minv_local; if(minv_local < minv)
minv = minv_local;
} }
if(maxv_local > maxv) if(maxv_local > maxv)
{ {
lock_guard<mutex> guard(max_mutex); lock_guard<mutex> guard(max_mutex);
maxv = maxv_local; if(maxv_local > maxv)
maxv = maxv_local;
} }
}); });
} }

7
libsrc/visualization/vssolution.hpp
View File

@ -43,7 +43,12 @@ class DLL_HEADER VisualSceneSolution : public VisualScene
Point<3> p; 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 surfellist;
int linelist; int linelist;
int element1dlist; int element1dlist;

31
ng/netgenpy.cpp
View File

@ -1,7 +1,7 @@
// a wrapper to load netgen-dll into python // a wrapper to load netgen-dll into python
#include <iostream> #include <iostream>
#include <boost/python.hpp> #include <../general/ngpython.hpp>
#ifdef WIN32 #ifdef WIN32
#define DLL_HEADER __declspec(dllimport) #define DLL_HEADER __declspec(dllimport)
@ -10,19 +10,26 @@
#endif #endif
void DLL_HEADER ExportNetgenMeshing(); void DLL_HEADER ExportNetgenMeshing(py::module &m);
void DLL_HEADER ExportMeshVis(); void DLL_HEADER ExportMeshVis(py::module &m);
void DLL_HEADER ExportCSG(); void DLL_HEADER ExportCSG(py::module &m);
void DLL_HEADER ExportCSGVis(); void DLL_HEADER ExportCSGVis(py::module &m);
void DLL_HEADER ExportGeom2d(); void DLL_HEADER ExportGeom2d(py::module &m);
BOOST_PYTHON_MODULE(libngpy) PYBIND11_PLUGIN(libngpy)
{ {
ExportCSG(); py::module ngpy("libngpy", "pybind netgen module");
ExportCSGVis(); py::module csg = ngpy.def_submodule("_csg", "pybind csg module");
ExportNetgenMeshing(); ExportCSG(csg);
ExportMeshVis(); py::module csgvis = ngpy.def_submodule("csgvis", "pybind csgvis module");
ExportGeom2d(); 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 // Force linking libnglib to libnetgenpy

8
ng/ngvisual.tcl
View File

@ -897,8 +897,8 @@ proc fieldlinesdialog { } {
set g $w.nb.settings set g $w.nb.settings
frame $g.linesettings -relief groove -borderwidth 3 ttk::frame $g.linesettings -relief groove -borderwidth 3
label $g.linesettings.title -text "\nLine Settings\n" ttk::label $g.linesettings.title -text "\nLine Settings\n"
# tixControl $g.linesettings.length -label "rel. Length: " -integer false \ # tixControl $g.linesettings.length -label "rel. Length: " -integer false \
# -variable visoptions.fieldlineslength -min 0.00001 -max 10000 -step 0.1 \ # -variable visoptions.fieldlineslength -min 0.00001 -max 10000 -step 0.1 \
# -options { # -options {
@ -951,8 +951,8 @@ proc fieldlinesdialog { } {
global visoptions.fieldlinestolerance global visoptions.fieldlinestolerance
frame $g.odesettings -relief groove -borderwidth 3 ttk::frame $g.odesettings -relief groove -borderwidth 3
label $g.odesettings.title -text "\nODE Settings\n" ttk::label $g.odesettings.title -text "\nODE Settings\n"
# tixControl $g.odesettings.tol -label "rel. Tolerance: " -integer false \ # tixControl $g.odesettings.tol -label "rel. Tolerance: " -integer false \
# -variable visoptions.fieldlinestolerance -min 0.00001 -max 1 -step 0.01 \ # -variable visoptions.fieldlinestolerance -min 0.00001 -max 1 -step 0.01 \
# -options { # -options {

76
ng/parallelfunc.cpp
View File

@ -138,7 +138,7 @@ void ParallelRun()
string redraw_cmd; string redraw_cmd;
// MyMPI_Recv (redraw_cmd, 0, MPI_TAG_VIS); // MyMPI_Recv (redraw_cmd, 0, MPI_TAG_VIS);
redraw_cmd = MyMPI_RecvCmd(); redraw_cmd = MyMPI_RecvCmd();
// PrintMessage (1, "Redraw - ", redraw_cmd); // PrintMessage (1, "Redraw - ", redraw_cmd);
static string displname; static string displname;
@ -147,17 +147,17 @@ void ParallelRun()
static Display * display = NULL; static Display * display = NULL;
static GLXContext context; static GLXContext context;
static XVisualInfo * visinfo = 0; static XVisualInfo * visinfo = 0;
// if (!display) // if (!display)
if (redraw_cmd == "init") if (redraw_cmd == "init")
{ {
MyMPI_Recv (displname, 0, MPI_TAG_VIS); MyMPI_Recv (displname, 0, MPI_TAG_VIS);
MyMPI_Recv (curDrawable, 0, MPI_TAG_VIS); MyMPI_Recv (curDrawable, 0, MPI_TAG_VIS);
MyMPI_Recv (contextid, 0, MPI_TAG_VIS); MyMPI_Recv (contextid, 0, MPI_TAG_VIS);
display = XOpenDisplay (displname.c_str()); display = XOpenDisplay (displname.c_str());
/* /*
PrintMessage (3, "displ - name = ", displname); PrintMessage (3, "displ - name = ", displname);
PrintMessage (3, "display = ", display, PrintMessage (3, "display = ", display,
@ -166,23 +166,27 @@ void ParallelRun()
" , h = ", XDisplayHeight (display, 0)); " , 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 cout << "P" << id << ": window-props: x = " << wx << ", y = " << wy
<< ", w = " << ww << ", h = " << wh << ", depth = " << depth << endl; << ", w = " << ww << ", h = " << wh << ", depth = " << depth << endl;
*/ */
#define VISUAL #define VISUAL
#ifdef VISUAL #ifdef VISUAL
#ifdef VISINFO_OLD
//this does not seem to work anymore (but might still be with togl1.7?)
// make a new GLXContext // make a new GLXContext
// first, generate a visual (copied from togl) // first, generate a visual (copied from togl)
@ -221,51 +225,68 @@ void ParallelRun()
attrib_list[attrib_count++] = GLX_DOUBLEBUFFER; attrib_list[attrib_count++] = GLX_DOUBLEBUFFER;
attrib_list[attrib_count++] = None; attrib_list[attrib_count++] = None;
visinfo = glXChooseVisual(display, 0, visinfo = glXChooseVisual(display, 0,
attrib_list); attrib_list);
cout << "have vis?" << endl;
if (visinfo) { if (visinfo) {
/* found a GLX visual! */ /* found a GLX visual! */
// cout << "found VISINFO !!!" << endl; // cout << "found VISINFO !!!" << endl;
cout << "found VISINFO !!!" << endl;
/* /*
int hi = 0; int hi = 0;
std::cout << "attribs = "; std::cout << "attribs = ";
while (attrib_list[hi] != None) while (attrib_list[hi] != None)
std::cout << attrib_list[hi++] << " "; std::cout << attrib_list[hi++] << " ";
std::cout << std::endl; std::cout << std::endl;
*/ */
break; break;
} }
} }
if (!visinfo) if (!visinfo)
cerr << "no VISINFO found" << endl; cerr << "no VISINFO found" << endl;
#else
//get all possible confs
int nconfs;
auto cptr = glXGetFBConfigs (display,0, &nconfs);
Array<int> conf_ids(nconfs);
for(int k=0;k<nconfs;k++)
glXGetFBConfigAttrib(display, cptr[k], GLX_FBCONFIG_ID, &(conf_ids[k]));
//get drawable->FBConfig->visual
unsigned int d_fbc_id;
glXQueryDrawable( display, curDrawable, GLX_FBCONFIG_ID, &d_fbc_id);
GLXFBConfig d_fbc;
for(int k=0;k<nconfs;k++)
if(d_fbc_id==conf_ids[k])
d_fbc = cptr[k];
visinfo = glXGetVisualFromFBConfig(display,d_fbc);
#endif
// context = glXCreateContext( display, visinfo, 0, /* curContext, */ False ); // context = glXCreateContext( display, visinfo, 0, /* curContext, */ False );
context = glXCreateContext( display, visinfo, glXImportContextEXT ( display, contextid ), False); context = glXCreateContext( display, visinfo, glXImportContextEXT ( display, contextid ), False);
// cout << "context = " << context << endl;
glXMakeCurrent (display, curDrawable, context); glXMakeCurrent (display, curDrawable, context);
#else #else
// try to get GLXcontext from the master. // try to get GLXcontext from the master.
// this needs an indirect context (BUT DOES NOT WORK ????) // this needs an indirect context (BUT DOES NOT WORK ????)
context = glXImportContextEXT ( display, contextid );
context = glXImportContextEXT ( display, contextid );
PrintMessage (1, "GLX-contextid = " , contextid, PrintMessage (1, "GLX-contextid = " , contextid,
" imported context ", context); " imported context ", context);
glXMakeCurrent (display, curDrawable, context); glXMakeCurrent (display, curDrawable, context);
#endif #endif
// PrintMessage (1, "redraw - init complete"); // PrintMessage (1, "redraw - init complete");
} }
if (redraw_cmd == "broadcast") if (redraw_cmd == "broadcast")
{ {
vsmesh.Broadcast (); vsmesh.Broadcast ();
@ -283,7 +304,6 @@ void ParallelRun()
vsmesh.BuildFilledList (false); vsmesh.BuildFilledList (false);
} }
if (redraw_cmd == "solsurfellist") if (redraw_cmd == "solsurfellist")
{ {
vssolution.DrawSurfaceElements(); vssolution.DrawSurfaceElements();

10
python/__init__.py
View File

@ -10,10 +10,14 @@ if __platform.startswith('darwin'):
# from libngpy import * # from libngpy import *
import libngpy import libngpy
# from libngpy import *
from . import csg # import libngpy
from . import meshing
# from . import csg
# from . import meshing
# from . import geom2d
del environ del environ
del path del path

11
python/csg.py
View File

@ -1,7 +1,8 @@
from libngpy.csg import * import libngpy
from libngpy._csg import *
import libngpy.csgvis as csgvis import libngpy.csgvis as csgvis
from libngpy.csgvis import MouseMove from libngpy.csgvis import MouseMove
from libngpy.meshing import MeshingParameters from libngpy._meshing import MeshingParameters
CSGeometry.VS = csgvis.VS CSGeometry.VS = csgvis.VS
@ -14,7 +15,11 @@ def VS (obj):
def csg_meshing_func (geom, **args): def csg_meshing_func (geom, **args):
return GenerateMesh (geom, MeshingParameters (**args)) if "mp" in args:
return GenerateMesh (geom, args["mp"])
else:
return GenerateMesh (geom, MeshingParameters (**args))
# return GenerateMesh (geom, MeshingParameters (**args))
CSGeometry.GenerateMesh = csg_meshing_func CSGeometry.GenerateMesh = csg_meshing_func

8
python/geom2d.py
View File

@ -1,5 +1,5 @@
from libngpy.geom2d import * from libngpy._geom2d import *
from libngpy.meshing import * from libngpy._meshing import *
tmp_generate_mesh = SplineGeometry.GenerateMesh tmp_generate_mesh = SplineGeometry.GenerateMesh
@ -42,8 +42,8 @@ def MakeCircle (geo, c, r, **args):
SplineGeometry.AddCircle = lambda geo, c, r, **args : MakeCircle(geo, c, r, **args) SplineGeometry.AddCircle = lambda geo, c, r, **args : MakeCircle(geo, c, r, **args)
SplineGeometry.AddRectangle = lambda geo, p1, p2, **args : MakeRectangle(geo, p1, p2, **args) SplineGeometry.AddRectangle = lambda geo, p1, p2, **args : MakeRectangle(geo, p1, p2, **args)
SplineGeometry.AddSegment = SplineGeometry.Append SplineGeometry.AddSegment = lambda *args, **kwargs : SplineGeometry.Append(*args, **kwargs)
SplineGeometry.AddPoint = SplineGeometry.AppendPoint SplineGeometry.AddPoint = lambda *args, **kwargs : SplineGeometry.AppendPoint(*args, **kwargs)
__all__ = ['SplineGeometry', 'unit_square'] __all__ = ['SplineGeometry', 'unit_square']

2
python/meshing.py
View File

@ -1 +1 @@
from libngpy.meshing import * from libngpy._meshing import *