diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 86105bfb..c5ee2e8c 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -14,8 +14,8 @@ stages: - x64 before_script: - "echo off" - - 'call "%VS140COMNTOOLS%\..\..\VC\bin\amd64\vcvars64.bat"' - - set CMAKE_GENERATOR=Visual Studio 14 2015 Win64 + - 'call "%VS2017INSTALLDIR%\VC\Auxiliary\Build\vcvars64"' + - set CMAKE_GENERATOR=Visual Studio 15 2017 Win64 - set CI_DIR=C:\ci\%CI_PIPELINE_ID% - set NETGEN_BUILD_DIR=%CI_DIR%\build - set INSTALL_DIR=%CI_DIR%\install @@ -100,6 +100,14 @@ test_guidelines: script: - docker run -e CCACHE_DIR=/ccache -v /mnt/ccache:/ccache netgen_${CI_BUILD_REF_NAME}:${UBUNTU_VERSION} bash /root/src/netgen/tests/build_guidelines.sh when: always + allow_failure: true + +# check if it compiles without spdlog +test_noSpdlog: + <<: *ubuntu + stage: test + script: + - docker run -e CCACHE_DIR=/ccache -v /mnt/ccache:/ccache netgen_${CI_BUILD_REF_NAME}:${UBUNTU_VERSION} bash /root/src/netgen/tests/build_nospdlog.sh cleanup_ubuntu: stage: cleanup diff --git a/CMakeLists.txt b/CMakeLists.txt index 14931c4d..bf7980e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,13 @@ if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE "RelWithDebInfo" CACHE STRING INTERNAL) endif(NOT CMAKE_BUILD_TYPE) -cmake_minimum_required(VERSION 3.1.3) +if(WIN32) + # we are linking to object libraries on Windows + cmake_minimum_required(VERSION 3.12) +else(WIN32) + cmake_minimum_required(VERSION 3.1.3) +endif(WIN32) + if(NOT WIN32) option( USE_NATIVE_ARCH "build which -march=native" ON) endif(NOT WIN32) @@ -19,6 +25,9 @@ option( USE_CCACHE "use ccache") option( USE_INTERNAL_TCL "Compile tcl files into the code and don't install them" ON) option( ENABLE_UNIT_TESTS "Enable Catch unit tests") option( ENABLE_CPP_CORE_GUIDELINES_CHECK "Enable cpp core guideline checks on ngcore" OFF) +option( USE_SPDLOG "Enable spd log logging" OFF) +option( DEBUG_LOG "Enable more debug output (may increase computation time) - only works with USE_SPDLOG=ON" OFF) +option( CHECK_RANGE "Check array range access, automatically enabled if built in debug mode" OFF) option( USE_SUPERBUILD "use ccache" ON) @@ -258,6 +267,7 @@ endif (USE_GUI) ####################################################################### if (USE_PYTHON) + add_subdirectory(external_dependencies/pybind11) add_definitions(-DNG_PYTHON) find_path(PYBIND_INCLUDE_DIR pybind11/pybind11.h HINTS ${PYTHON_INCLUDE_DIR}) if( PYBIND_INCLUDE_DIR ) @@ -350,6 +360,11 @@ endif(ENABLE_UNIT_TESTS) ####################################################################### +if(USE_SPDLOG) + include(${CMAKE_CURRENT_SOURCE_DIR}/cmake/external_projects/spdlog.cmake) + include_directories(${SPDLOG_INCLUDE_DIR}) +endif(USE_SPDLOG) + if(ENABLE_CPP_CORE_GUIDELINES_CHECK) find_program( CLANG_TIDY_EXE diff --git a/cmake/SuperBuild.cmake b/cmake/SuperBuild.cmake index 36ea3fcb..646de088 100644 --- a/cmake/SuperBuild.cmake +++ b/cmake/SuperBuild.cmake @@ -142,6 +142,9 @@ set_vars( NETGEN_CMAKE_ARGS CMAKE_INSTALL_PREFIX ENABLE_UNIT_TESTS ENABLE_CPP_CORE_GUIDELINES_CHECK + USE_SPDLOG + DEBUG_LOG + CHECK_RANGE ) # propagate all variables set on the command line using cmake -DFOO=BAR diff --git a/cmake/external_projects/spdlog.cmake b/cmake/external_projects/spdlog.cmake new file mode 100644 index 00000000..d932ebef --- /dev/null +++ b/cmake/external_projects/spdlog.cmake @@ -0,0 +1,18 @@ +include(ExternalProject) +find_program(GIT_EXECUTABLE git) + +ExternalProject_Add( + project_spdlog + PREFIX ${CMAKE_BINARY_DIR}/spdlog + GIT_REPOSITORY https://github.com/gabime/spdlog.git + GIT_TAG v1.2.1 + TIMEOUT 01 + UPDATE_COMMAND "" + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + LOG_DOWNLOAD ON + ) + +ExternalProject_Get_Property(project_spdlog source_dir) +set(SPDLOG_INCLUDE_DIR ${source_dir}/include) diff --git a/cmake/external_projects/tcltk.cmake b/cmake/external_projects/tcltk.cmake index b7f0ba85..b071e495 100644 --- a/cmake/external_projects/tcltk.cmake +++ b/cmake/external_projects/tcltk.cmake @@ -1,6 +1,6 @@ if(APPLE) # use system tcl/tk - if(${PYTHON_VERSION_STRING} STREQUAL "3.7") + if((${PYTHON_VERSION_STRING} VERSION_EQUAL "3.7") OR (${PYTHON_VERSION_STRING} VERSION_GREATER "3.7")) # fetch tcl/tk sources to match the one used in Python 3.7 ExternalProject_Add(project_tcl URL "https://prdownloads.sourceforge.net/tcl/tcl8.6.8-src.tar.gz" diff --git a/libsrc/core/.clang-tidy b/libsrc/core/.clang-tidy index 290188fb..5742a3d8 100644 --- a/libsrc/core/.clang-tidy +++ b/libsrc/core/.clang-tidy @@ -1,5 +1,5 @@ -Checks: '*,-clang-analyzer-alpha.*,-*braces-around-statements,-fuchsia-*,-google-runtime-references,-readability-implicit-bool-conversion,-google-explicit-constructor,-hicpp-explicit-conversions,-google-runtime-int,-llvm-header-guard' +Checks: '*,-clang-analyzer-alpha.*,-*braces-around-statements,-fuchsia-*,-google-runtime-references,-readability-implicit-bool-conversion,-google-explicit-constructor,-hicpp-explicit-conversions,-google-runtime-int,-llvm-header-guard,-modernize-pass-by-value' CheckOptions: - key: cppcoreguidelines-special-member-functions.AllowSoleDefaultDtor value: 1 -WarningsAsErrors: '*' \ No newline at end of file +WarningsAsErrors: '*' diff --git a/libsrc/core/CMakeLists.txt b/libsrc/core/CMakeLists.txt index 92f2e8b7..46a2230b 100644 --- a/libsrc/core/CMakeLists.txt +++ b/libsrc/core/CMakeLists.txt @@ -1,16 +1,49 @@ -add_library(ngcore SHARED archive.cpp) -target_compile_definitions(ngcore PRIVATE -DNGCORE_EXPORTS) +add_library(ngcore SHARED archive.cpp logging.cpp paje_trace.cpp utils.cpp profiler.cpp) + +target_compile_definitions(ngcore PRIVATE NGCORE_EXPORTS) +if(NOT WIN32) + target_compile_options(ngcore PRIVATE -fvisibility=hidden) +endif(NOT WIN32) + +target_compile_definitions(ngcore PUBLIC $<$:NETGEN_ENABLE_CHECK_RANGE>) + +if(CHECK_RANGE) + target_compile_definitions(ngcore PUBLIC NETGEN_ENABLE_CHECK_RANGE) +endif(CHECK_RANGE) + +if(USE_SPDLOG) + include_directories(${SPDLOG_INCLUDE_DIR}) + install(DIRECTORY ${SPDLOG_INCLUDE_DIR} + DESTINATION ${NG_INSTALL_DIR_INCLUDE} + ) + add_dependencies(ngcore project_spdlog) + target_compile_definitions(ngcore PUBLIC NETGEN_USE_SPDLOG) + if(DEBUG_LOG) + target_compile_definitions(ngcore PUBLIC NETGEN_LOG_DEBUG) + endif(DEBUG_LOG) +endif(USE_SPDLOG) install(TARGETS ngcore DESTINATION ${NG_INSTALL_DIR} COMPONENT netgen) if(USE_PYTHON) + target_compile_definitions(ngcore PUBLIC NETGEN_PYTHON) target_include_directories(ngcore PUBLIC ${PYTHON_INCLUDE_DIRS}) + target_link_libraries(ngcore PUBLIC ${PYTHON_LIBRARIES}) endif(USE_PYTHON) -install(FILES ngcore.hpp archive.hpp type_traits.hpp version.hpp ngcore_api.hpp +install(FILES ngcore.hpp archive.hpp type_traits.hpp version.hpp ngcore_api.hpp logging.hpp + exception.hpp symboltable.hpp paje_trace.hpp utils.hpp profiler.hpp DESTINATION ${NG_INSTALL_DIR_INCLUDE}/core COMPONENT netgen_devel) if(ENABLE_CPP_CORE_GUIDELINES_CHECK) set_target_properties(ngcore PROPERTIES CXX_CLANG_TIDY "${DO_CLANG_TIDY}") endif(ENABLE_CPP_CORE_GUIDELINES_CHECK) + +if(USE_PYTHON) + pybind11_add_module(pyngcore SHARED python_ngcore.cpp) + target_link_libraries(pyngcore PUBLIC ngcore ${PYTHON_LIBRARIES}) + set_target_properties(pyngcore PROPERTIES INSTALL_RPATH "${NG_RPATH_TOKEN}/${NETGEN_PYTHON_RPATH}") + install(TARGETS pyngcore DESTINATION ${NG_INSTALL_DIR_PYTHON} COMPONENT netgen) +endif(USE_PYTHON) + diff --git a/libsrc/core/archive.cpp b/libsrc/core/archive.cpp index 64789904..9d7f2cd5 100644 --- a/libsrc/core/archive.cpp +++ b/libsrc/core/archive.cpp @@ -19,16 +19,6 @@ namespace ngcore void SetLibraryVersion(const std::string& library, const VersionInfo& version) { library_versions[library] = version; } -#ifdef WIN32 - // windows does demangling in typeid(T).name() - std::string Demangle(const char* typeinfo) { return typeinfo; } -#else - std::string Demangle(const char* typeinfo) { int status; return abi::__cxa_demangle(typeinfo, - nullptr, - nullptr, - &status); } -#endif - // clang-tidy should ignore this static object static std::unique_ptr> type_register; // NOLINT const detail::ClassArchiveInfo& Archive :: GetArchiveRegister(const std::string& classname) diff --git a/libsrc/core/archive.hpp b/libsrc/core/archive.hpp index fd2ee275..c1182b40 100644 --- a/libsrc/core/archive.hpp +++ b/libsrc/core/archive.hpp @@ -3,31 +3,32 @@ #include // for complex #include // for size_t, strlen -#include // for operator<<, ifstream, ofstream, basic... +#include // for ifstream, ofstream #include // for function -#include // for map, _Rb_tree_iterator -#include // for __shared_ptr_access, __shared_ptr_acc... -#include // for runtime_error -#include // for string, operator+ +#include // for map +#include // for shared_ptr +#include // for string #include // for declval, enable_if, false_type, is_co... #include // for type_info #include // for move, swap, pair #include // for vector -#include "ngcore_api.hpp" // for NGCORE_API, unlikely +#include "exception.hpp" // for UnreachableCodeException, Exception +#include "logging.hpp" // for logger +#include "ngcore_api.hpp" // for NGCORE_API #include "type_traits.hpp" // for all_of_tmpl +#include "utils.hpp" // for Demangle, unlikely #include "version.hpp" // for VersionInfo -#ifdef NG_PYTHON +#ifdef NETGEN_PYTHON #include -#endif // NG_PYTHON +#endif // NETGEN_PYTHON namespace ngcore { // Libraries using this archive can store their version here to implement backwards compatibility NGCORE_API const VersionInfo& GetLibraryVersion(const std::string& library); NGCORE_API void SetLibraryVersion(const std::string& library, const VersionInfo& version); - NGCORE_API std::string Demangle(const char* typeinfo); class NGCORE_API Archive; @@ -36,7 +37,7 @@ namespace ngcore // create new pointer of type T if it is default constructible, else throw template T* constructIfPossible_impl(Rest... /*unused*/) - { throw std::runtime_error(std::string(Demangle(typeid(T).name())) + " is not default constructible!"); } + { throw Exception(std::string(Demangle(typeid(T).name())) + " is not default constructible!"); } template::value>::type> T* constructIfPossible_impl(int /*unused*/) { return new T; } // NOLINT @@ -96,28 +97,34 @@ namespace ngcore { const bool is_output; // how many different shared_ptr/pointer have been (un)archived - int shared_ptr_count, ptr_count; + int shared_ptr_count{0}, ptr_count{0}; // maps for archived shared pointers and pointers - std::map shared_ptr2nr, ptr2nr; + std::map shared_ptr2nr{}, ptr2nr{}; // vectors for storing the unarchived (shared) pointers - std::vector> nr2shared_ptr; - std::vector nr2ptr; + std::vector> nr2shared_ptr{}; + std::vector nr2ptr{}; protected: bool shallow_to_python = false; + std::map version_map = GetLibraryVersions(); + std::shared_ptr logger = GetLogger("Archive"); public: Archive() = delete; Archive(const Archive&) = delete; Archive(Archive&&) = delete; - Archive (bool ais_output) : - is_output(ais_output), shared_ptr_count(0), ptr_count(0) { ; } + Archive (bool ais_output) : is_output(ais_output) { ; } virtual ~Archive() { ; } + // If the object is pickled, all shallow archived objects will be pickled as a list, + // instead of written as a binary archive. This allows pickle to serialize every object only + // once and put them together correctly afterwards. Therefore all objects that may live in + // Python should be archived using this Shallow function. If Shallow is called from C++ code + // it archives the object normally. template Archive& Shallow(T& val) { static_assert(detail::is_any_pointer, "ShallowArchive must be given pointer type!"); -#ifdef NG_PYTHON +#ifdef NETGEN_PYTHON if(shallow_to_python) { if(is_output) @@ -126,25 +133,28 @@ namespace ngcore val = pybind11::cast(ShallowInPython()); } else -#endif // NG_PYTHON +#endif // NETGEN_PYTHON *this & val; return *this; } -#ifdef NG_PYTHON - virtual void ShallowOutPython(pybind11::object /*unused*/) // NOLINT (copy by val is ok for this virt func) - { throw std::runtime_error("Should not get in ShallowToPython base class implementation!"); } +#ifdef NETGEN_PYTHON + virtual void ShallowOutPython(const pybind11::object& /*unused*/) + { throw UnreachableCodeException{}; } virtual pybind11::object ShallowInPython() - { throw std::runtime_error("Should not get in ShallowFromPython base class implementation!"); } -#endif // NG_PYTHON + { throw UnreachableCodeException{}; } +#endif // NETGEN_PYTHON Archive& operator=(const Archive&) = delete; Archive& operator=(Archive&&) = delete; bool Output () const { return is_output; } bool Input () const { return !is_output; } - virtual const VersionInfo& GetVersion(const std::string& library) - { return GetLibraryVersions()[library]; } + const VersionInfo& GetVersion(const std::string& library) + { return version_map[library]; } + + // only used for PyArchive + virtual void NeedsVersion(const std::string& /*unused*/, const std::string& /*unused*/) {} // Pure virtual functions that have to be implemented by In-/OutArchive virtual Archive & operator & (double & d) = 0; @@ -157,7 +167,7 @@ namespace ngcore virtual Archive & operator & (std::string & str) = 0; virtual Archive & operator & (char *& str) = 0; - virtual Archive & operator & (VersionInfo & version) + Archive & operator & (VersionInfo & version) { if(Output()) (*this) << version.to_string(); @@ -198,6 +208,47 @@ namespace ngcore Do(&v[0], size); return (*this); } + + // archive implementation for enums + template + auto operator & (T& val) -> typename std::enable_if::value, Archive&>::type + { + int enumval; + if(Output()) + enumval = int(val); + *this & enumval; + if(Input()) + val = T(enumval); + return *this; + } + + // vector has special implementation (like a bitarray) therefore + // it needs a special overload (this could probably be more efficient, but we + // don't use it that often anyway) + Archive& operator& (std::vector& v) + { + logger->debug("In special archive for std::vector"); + size_t size; + if(Output()) + size = v.size(); + (*this) & size; + if(Input()) + { + v.resize(size); + bool b; + for(size_t i=0; i Archive& operator& (std::map& map) { @@ -261,28 +312,40 @@ namespace ngcore { if(Output()) { + logger->debug("Store shared ptr of type {}", Demangle(typeid(T).name())); // save -2 for nullptr if(!ptr) - return (*this) << -2; + { + logger->debug("Storing nullptr"); + return (*this) << -2; + } void* reg_ptr = ptr.get(); bool neededDowncast = false; // Downcasting is only possible for our registered classes if(typeid(T) != typeid(*ptr)) { + logger->debug("Typids are different: {} vs {}", + Demangle(typeid(T).name()), + Demangle(typeid(*ptr).name())); if(!IsRegistered(Demangle(typeid(*ptr).name()))) - throw std::runtime_error(std::string("Archive error: Polymorphic type ") - + Demangle(typeid(*ptr).name()) - + " not registered for archive"); + throw Exception(std::string("Archive error: Polymorphic type ") + + Demangle(typeid(*ptr).name()) + + " not registered for archive"); reg_ptr = GetArchiveRegister(Demangle(typeid(*ptr).name())).downcaster(typeid(T), ptr.get()); // if there was a true downcast we have to store more information - if(reg_ptr != static_cast(ptr.get()) ) + if(reg_ptr != static_cast(ptr.get())) + { + logger->debug("Multiple/Virtual inheritance involved, need to cast pointer"); neededDowncast = true; + } } auto pos = shared_ptr2nr.find(reg_ptr); // if not found store -1 and the pointer if(pos == shared_ptr2nr.end()) { + logger->debug("Didn't find the shared_ptr, create new registry entry at {}", + shared_ptr_count); auto p = ptr.get(); (*this) << -1; (*this) & neededDowncast & p; @@ -293,23 +356,27 @@ namespace ngcore return *this; } // if found store the position and if it has to be downcasted and how + logger->debug("Found shared_ptr at position {}", pos->second); (*this) << pos->second << neededDowncast; if(neededDowncast) (*this) << Demangle(typeid(*ptr).name()); } else // Input { + logger->debug("Reading shared_ptr of type {}", Demangle(typeid(T).name())); int nr; (*this) & nr; // -2 restores a nullptr if(nr == -2) { + logger->debug("Reading a nullptr"); ptr = nullptr; return *this; } // -1 restores a new shared ptr by restoring the inner pointer and creating a shared_ptr to it if (nr == -1) { + logger->debug("Createing new shared_ptr"); T* p = nullptr; bool neededDowncast; (*this) & neededDowncast & p; @@ -317,6 +384,7 @@ namespace ngcore // if we did downcast we need to store a shared_ptr to the true object if(neededDowncast) { + logger->debug("Shared pointer needed downcasting"); std::string name; (*this) & name; auto info = GetArchiveRegister(name); @@ -327,15 +395,20 @@ namespace ngcore ptr.get()))); } else - nr2shared_ptr.push_back(ptr); + { + logger->debug("Shared pointer didn't need downcasting"); + nr2shared_ptr.push_back(ptr); + } } else { + logger->debug("Reading already existing pointer at entry {}", nr); auto other = nr2shared_ptr[nr]; bool neededDowncast; (*this) & neededDowncast; if(neededDowncast) { + logger->debug("Shared pointer needed pointer downcast"); // if there was a downcast we can expect the class to be registered (since archiving // wouldn't have worked else) std::string name; @@ -348,7 +421,10 @@ namespace ngcore other.get()))); } else - ptr = std::static_pointer_cast(other); + { + logger->debug("Shared pointer didn't need pointer casts"); + ptr = std::static_pointer_cast(other); + } } } return *this; @@ -360,31 +436,45 @@ namespace ngcore { if (Output()) { + logger->debug("Store pointer of type {}",Demangle(typeid(T).name())); // if the pointer is null store -2 if (!p) - return (*this) << -2; + { + logger->debug("Storing nullptr"); + return (*this) << -2; + } auto reg_ptr = static_cast(p); if(typeid(T) != typeid(*p)) { + logger->debug("Typeids are different: {} vs {}", + Demangle(typeid(T).name()), + Demangle(typeid(*p).name())); if(!IsRegistered(Demangle(typeid(*p).name()))) - throw std::runtime_error(std::string("Archive error: Polymorphic type ") - + Demangle(typeid(*p).name()) - + " not registered for archive"); + throw Exception(std::string("Archive error: Polymorphic type ") + + Demangle(typeid(*p).name()) + + " not registered for archive"); reg_ptr = GetArchiveRegister(Demangle(typeid(*p).name())).downcaster(typeid(T), static_cast(p)); + if(reg_ptr != static_cast(p)) + { + logger->debug("Multiple/Virtual inheritance involved, need to cast pointer"); + } } auto pos = ptr2nr.find(reg_ptr); // if the pointer is not found in the map create a new entry if (pos == ptr2nr.end()) { + logger->debug("Didn't find pointer, create new registry entry at {}", + ptr_count); ptr2nr[reg_ptr] = ptr_count++; if(typeid(*p) == typeid(T)) if (std::is_constructible::value) { + logger->debug("Store standard class pointer (no virt. inh,...)"); return (*this) << -1 & (*p); } else - throw std::runtime_error(std::string("Archive error: Class ") + - Demangle(typeid(*p).name()) + " does not provide a default constructor!"); + throw Exception(std::string("Archive error: Class ") + + Demangle(typeid(*p).name()) + " does not provide a default constructor!"); else { // if a pointer to a base class is archived, the class hierarchy must be registered @@ -392,9 +482,10 @@ namespace ngcore // implement a void DoArchive(Archive&) member function // To recreate the object we need to store the true type of it if(!IsRegistered(Demangle(typeid(*p).name()))) - throw std::runtime_error(std::string("Archive error: Polymorphic type ") - + Demangle(typeid(*p).name()) - + " not registered for archive"); + throw Exception(std::string("Archive error: Polymorphic type ") + + Demangle(typeid(*p).name()) + + " not registered for archive"); + logger->debug("Store a possibly more complicated pointer"); return (*this) << -3 << Demangle(typeid(*p).name()) & (*p); } } @@ -402,27 +493,37 @@ namespace ngcore { (*this) & pos->second; bool downcasted = !(reg_ptr == static_cast(p) ); + logger->debug("Store a the existing position in registry at {}", pos->second); + logger->debug("Pointer {} downcasting", downcasted ? "needs" : "doesn't need"); // store if the class has been downcasted and the name (*this) << downcasted << Demangle(typeid(*p).name()); } } else { + logger->debug("Reading pointer of type {}", Demangle(typeid(T).name())); int nr; (*this) & nr; if (nr == -2) // restore a nullptr + { + logger->debug("Loading a nullptr"); p = nullptr; + } else if (nr == -1) // create a new pointer of standard type (no virtual or multiple inheritance,...) { + logger->debug("Load a new pointer to a simple class"); p = detail::constructIfPossible(); nr2ptr.push_back(p); (*this) & *p; } else if(nr == -3) // restore one of our registered classes that can have multiple inheritance,... { + logger->debug("Load a new pointer to a potentially more complicated class " + "(allows for multiple/virtual inheritance,...)"); // As stated above, we want this special behaviour only for our classes that implement DoArchive std::string name; (*this) & name; + logger->debug("Name = {}", name); auto info = GetArchiveRegister(name); // the creator creates a new object of type name, and returns a void* pointing // to T (which may have an offset) @@ -434,9 +535,11 @@ namespace ngcore } else { + logger->debug("Restoring pointer to already existing object at registry position {}", nr); bool downcasted; std::string name; (*this) & downcasted & name; + logger->debug("{} object of type {}", downcasted ? "Downcasted" : "Not downcasted", name); if(downcasted) { // if the class has been downcasted we can assume it is in the register @@ -491,11 +594,11 @@ namespace ngcore { static void* tryUpcast (const std::type_info& /*unused*/, T* /*unused*/) { - throw std::runtime_error("Upcast not successful, some classes are not registered properly for archiving!"); + throw Exception("Upcast not successful, some classes are not registered properly for archiving!"); } static void* tryDowncast (const std::type_info& /*unused*/, void* /*unused*/) { - throw std::runtime_error("Downcast not successful, some classes are not registered properly for archiving!"); + throw Exception("Downcast not successful, some classes are not registered properly for archiving!"); } }; @@ -507,7 +610,7 @@ namespace ngcore try { return GetArchiveRegister(Demangle(typeid(B1).name())). upcaster(ti, static_cast(dynamic_cast(p))); } - catch(std::exception&) + catch(const Exception&) { return Caster::tryUpcast(ti, p); } } @@ -520,7 +623,7 @@ namespace ngcore return dynamic_cast(static_cast(GetArchiveRegister(Demangle(typeid(B1).name())). downcaster(ti, p))); } - catch(std::exception&) + catch(const Exception&) { return Caster::tryDowncast(ti, p); } @@ -536,7 +639,7 @@ namespace ngcore { static_assert(detail::all_of_tmpl::value...>, "Variadic template arguments must be base classes of T"); - detail::ClassArchiveInfo info; + detail::ClassArchiveInfo info {}; info.creator = [this,&info](const std::type_info& ti) -> void* { return typeid(T) == ti ? detail::constructIfPossible() : Archive::Caster::tryUpcast(ti, detail::constructIfPossible()); }; @@ -805,7 +908,7 @@ namespace ngcore } }; -#ifdef NG_PYTHON +#ifdef NETGEN_PYTHON template class PyArchive : public ARCHIVE @@ -813,7 +916,12 @@ namespace ngcore private: pybind11::list lst; size_t index = 0; + std::map version_needed; + protected: using ARCHIVE::stream; + using ARCHIVE::version_map; + using ARCHIVE::logger; + using ARCHIVE::GetLibraryVersions; public: PyArchive(const pybind11::object& alst = pybind11::none()) : ARCHIVE(std::make_shared()), @@ -821,8 +929,30 @@ namespace ngcore { ARCHIVE::shallow_to_python = true; if(Input()) - stream = std::make_shared - (pybind11::cast(lst[pybind11::len(lst)-1])); + { + stream = std::make_shared + (pybind11::cast(lst[pybind11::len(lst)-1])); + *this & version_needed; + logger->debug("versions needed for unpickling = {}", version_needed); + for(auto& libversion : version_needed) + if(libversion.second > GetLibraryVersion(libversion.first)) + throw Exception("Error in unpickling data:\nLibrary " + libversion.first + + " must be at least " + libversion.second.to_string()); + stream = std::make_shared + (pybind11::cast(lst[pybind11::len(lst)-2])); + *this & version_map; + stream = std::make_shared + (pybind11::cast(lst[pybind11::len(lst)-3])); + } + } + + void NeedsVersion(const std::string& library, const std::string& version) override + { + if(Output()) + { + logger->debug("Need version {} of library {}.", version, library); + version_needed[library] = version_needed[library] > version ? version_needed[library] : version; + } } using ARCHIVE::Output; @@ -831,11 +961,20 @@ namespace ngcore using ARCHIVE::operator&; using ARCHIVE::operator<<; using ARCHIVE::GetVersion; - void ShallowOutPython(pybind11::object val) override { lst.append(val); } + void ShallowOutPython(const pybind11::object& val) override { lst.append(val); } pybind11::object ShallowInPython() override { return lst[index++]; } pybind11::list WriteOut() { + FlushBuffer(); + lst.append(pybind11::bytes(std::static_pointer_cast(stream)->str())); + stream = std::make_shared(); + *this & GetLibraryVersions(); + FlushBuffer(); + lst.append(pybind11::bytes(std::static_pointer_cast(stream)->str())); + stream = std::make_shared(); + logger->debug("Writeout version needed = {}", version_needed); + *this & version_needed; FlushBuffer(); lst.append(pybind11::bytes(std::static_pointer_cast(stream)->str())); return lst; @@ -843,27 +982,31 @@ namespace ngcore }; template - auto NGSPickle(bool printoutput=false) + auto NGSPickle() { - return pybind11::pickle([printoutput](T* self) + return pybind11::pickle([](T* self) { PyArchive ar; ar & self; auto output = pybind11::make_tuple(ar.WriteOut()); - if(printoutput) - pybind11::print("pickle output of", Demangle(typeid(T).name()),"=", output); + GetLogger("Archive")->trace("Pickling output for object of type {} = {}", + Demangle(typeid(T).name()), + std::string(pybind11::str(output))); return output; }, [](pybind11::tuple state) { T* val = nullptr; + GetLogger("Archive")->trace("State for unpickling of object of type {} = {}", + Demangle(typeid(T).name()), + std::string(pybind11::str(state[0]))); PyArchive ar(state[0]); ar & val; return val; }); } -#endif // NG_PYTHON +#endif // NETGEN_PYTHON } // namespace ngcore #endif // NETGEN_CORE_ARCHIVE_HPP diff --git a/libsrc/core/exception.hpp b/libsrc/core/exception.hpp new file mode 100644 index 00000000..35359e65 --- /dev/null +++ b/libsrc/core/exception.hpp @@ -0,0 +1,90 @@ +#ifndef NETGEN_CORE_EXCEPTION_HPP +#define NETGEN_CORE_EXCEPTION_HPP + +#include // for stringstream +#include // for exception +#include // for string + +#include "ngcore_api.hpp" // for NGCORE_API + +namespace ngcore +{ + // Exception for code that shouldn't be executed + class NGCORE_API UnreachableCodeException : public std::exception + { + const char* what() const noexcept override + { + return "Shouldn't get here, something went wrong!"; + } + }; + + // Default exception class + class NGCORE_API Exception : public std::exception + { + /// a verbal description of the exception + std::string m_what; + public: + Exception() = default; + Exception(const Exception&) = default; + Exception(Exception&&) = default; + Exception(const std::string& s) : m_what(s) {} + Exception(const char* s) : m_what(s) {} + ~Exception() override = default; + + Exception& operator =(const Exception&) = default; + Exception& operator =(Exception&&) noexcept = default; + + /// append string to description + Exception & Append (const std::string & s) { m_what += s; return *this; } + /// append string to description + Exception & Append (const char * s) { m_what += s; return *this; } + + /// verbal description of exception + const std::string & What() const { return m_what; } + + /// implement virtual function of std::exception + const char* what() const noexcept override { return m_what.c_str(); } + }; + + // Out of Range exception + class NGCORE_API RangeException : public Exception + { + public: + /// where it occurs, index, minimal and maximal indices + RangeException (const std::string & where, + int ind, int imin, int imax) : Exception("") + { + std::stringstream str; + str << where << ": index " << ind << " out of range [" << imin << "," << imax << "]\n"; + Append (str.str()); + } + + template + RangeException(const std::string& where, const T& value) + { + std::stringstream str; + str << where << " called with wrong value " << value << "\n"; + Append(str.str()); + } + }; + + // Exception used if no simd implementation is available to fall back to standard evaluation + class NGCORE_API ExceptionNOSIMD : public Exception + { public: using Exception::Exception; }; +} // namespace ngcore + +#define NETGEN_CORE_NGEXEPTION_STR_HELPER(x) #x +#define NETGEN_CORE_NGEXEPTION_STR(x) NETGEN_CORE_NGEXEPTION_STR_HELPER(x) + +// Convenience macro to append file name and line of exception origin to the string +#define NG_EXCEPTION(s) ngcore::Exception(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t"+std::string(s)) + +#ifdef NETGEN_ENABLE_CHECK_RANGE +#define NETGEN_CHECK_RANGE(value, min, max) \ + { if ((value)<(min) || (value)>=(max)) \ + throw ngcore::RangeException(__FILE__ ":" NETGEN_CORE_NGEXEPTION_STR(__LINE__) "\t", (value), (min), (max)); } +#else // NETGEN_ENABLE_CHECK_RANGE +#define NETGEN_CHECK_RANGE(value, min, max) +#endif // NETGEN_ENABLE_CHECK_RANGE + +#endif // NETGEN_CORE_EXCEPTION_HPP diff --git a/libsrc/core/logging.cpp b/libsrc/core/logging.cpp new file mode 100644 index 00000000..eaaedf02 --- /dev/null +++ b/libsrc/core/logging.cpp @@ -0,0 +1,137 @@ +#include "logging.hpp" + +#ifdef NETGEN_USE_SPDLOG + +#include +#include +#include + +#else // NETGEN_USE_SPDLOG +#include +#endif // NETGEN_USE_SPDLOG + + +namespace ngcore +{ + + void Logger::log(level::level_enum level, std::string && s) + { +#ifdef NETGEN_USE_SPDLOG + logger->log(spdlog::level::level_enum(level), s); +#else // NETGEN_USE_SPDLOG + if(level>level::debug) + std::clog << s << '\n'; +#endif // NETGEN_USE_SPDLOG + } + +#ifdef NETGEN_USE_SPDLOG + namespace detail + { + std::vector>& GetDefaultSinks() + { + static std::vector> sinks = + { std::make_shared() }; + return sinks; + } + std::shared_ptr CreateDefaultLogger(const std::string& name) + { + auto& default_sinks = GetDefaultSinks(); + auto logger = std::make_shared(name, default_sinks.begin(), default_sinks.end()); + spdlog::details::registry::instance().register_and_init(logger); + return logger; + } + } // namespace detail + + std::shared_ptr GetLogger(const std::string& name) + { + auto logger = spdlog::get(name); + if(!logger) + logger = detail::CreateDefaultLogger(name); + return std::make_shared(logger); + } + + void SetLoggingLevel(spdlog::level::level_enum level, const std::string& name) + { + if(!name.empty()) + spdlog::get(name)->set_level(level); + else + spdlog::set_level(level); + } + + void AddFileSink(const std::string& filename, spdlog::level::level_enum level, const std::string& logger) + { + auto sink = std::make_shared(filename); + sink->set_level(level); + if(!logger.empty()) + GetLogger(logger)->logger->sinks().push_back(sink); + else + { + detail::GetDefaultSinks().push_back(sink); + spdlog::details::registry::instance().apply_all([sink](auto logger) { logger->sinks().push_back(sink); }); + } + } + + void AddConsoleSink(spdlog::level::level_enum level, const std::string& logger) + { + auto sink = std::make_shared(); + sink->set_level(level); + if(!logger.empty()) + GetLogger(logger)->logger->sinks().push_back(sink); + else + { + detail::GetDefaultSinks().push_back(sink); + spdlog::details::registry::instance().apply_all([sink](auto logger) { logger->sinks().push_back(sink); }); + } + } + + void ClearLoggingSinks(const std::string& logger) + { + if(!logger.empty()) + GetLogger(logger)->logger->sinks().clear(); + else + { + detail::GetDefaultSinks().clear(); + spdlog::details::registry::instance().apply_all([](auto logger) { logger->sinks().clear(); }); + } + } + + void FlushOnLoggingLevel(spdlog::level::level_enum level, const std::string& logger) + { + if(!logger.empty()) + GetLogger(logger)->logger->flush_on(level); + else + spdlog::flush_on(level); + } + +#else // NETGEN_USE_SPDLOG +} //namespace ngcore + +namespace spdlog +{ + class logger + { + public: + logger() = default; + }; +} // namespace spdlog + +namespace ngcore +{ + + // Dummy functions if no spdlog is available + + std::shared_ptr GetLogger(const std::string& /*unused*/) + { + return std::make_shared(std::make_shared()); + } + + void SetLoggingLevel(level::level_enum /*unused*/, const std::string& /*unused*/) {} + void AddFileSink(const std::string& /*unused*/, level::level_enum /*unused*/, + const std::string& /*unused*/) + {} + void AddConsoleSink(level::level_enum /*unused*/, const std::string& /*unused*/) {} + void ClearLoggingSinks(const std::string& /*unused*/) {} + void FlushOnLoggingLevel(level::level_enum /*unused*/, const std::string& /*unused*/) {} +} //namespace ngcore + +#endif // NETGEN_USE_SPDLOG diff --git a/libsrc/core/logging.hpp b/libsrc/core/logging.hpp new file mode 100644 index 00000000..1553ad54 --- /dev/null +++ b/libsrc/core/logging.hpp @@ -0,0 +1,124 @@ +#ifndef NETGEN_CORE_LOGGING_HPP +#define NETGEN_CORE_LOGGING_HPP + +#undef NETGEN_USE_SPDLOG +#include +#include +#include +#include + +#include "exception.hpp" +#include "ngcore_api.hpp" +#include "utils.hpp" + +#ifdef NETGEN_USE_SPDLOG +#include +#include // to be able to parse anything to logger that implements operator << ostream +#ifdef NETGEN_LOG_DEBUG +#define SPDLOG_DEBUG_ON +#define NETGEN_DEBUG_LOG(logger, ...) SPDLOG_DEBUG(logger, __VA_ARGS__) +#endif // NETGEN_LOG_DEBUG +#endif // NETGEN_USE_SPDLOG + +#ifndef NETGEN_DEBUG_LOG +#define NETGEN_DEBUG_LOG(logger, ...) +#endif // NETGEN_DEBUG_LOG + +namespace spdlog +{ + class logger; +} // namespace spdlog + +namespace ngcore +{ + namespace level + { + enum level_enum + { + trace = 0, + debug = 1, + info = 2, + warn = 3, + err = 4, + critical = 5, + off = 6 + }; + } // namespace level + + class Logger + { + public: + std::shared_ptr logger; + + Logger(std::shared_ptr l) : logger(std::move(l)) {} + + void NGCORE_API log( level::level_enum level, std::string && s); + +#ifdef NETGEN_USE_SPDLOG + template + void log( level::level_enum level, const char* str, Args ... args) + { + log(level, fmt::format(str, args...)); + } +#else // NETGEN_USE_SPDLOG + template + std::string replace(std::string s, const T & t) + { + auto p0 = s.find_first_of('{'); + auto p1 = s.find_first_of('}', p0); + if(p0==std::string::npos || p1==std::string::npos) + throw Exception("invalid format string"); + s.replace(p0, p1-p0+1, ToString(t)); + return s; + } + + std::string log_helper(std::string s) + { + return s; + } + + template + std::string log_helper(std::string s, const T &t) + { + return replace(s,t); + } + + template + std::string log_helper( std::string s, const T &t, Args ... args) + { + return log_helper(replace(s,t), args...); + } + + template + void log( level::level_enum level, const char* str, Args ... args) + { + log(level, log_helper(std::string(str), args...)); + } +#endif // NETGEN_USE_SPDLOG + + template + void trace( const char* str, Args ... args) { log(level::level_enum::trace, str, args...); } + template + void debug( const char* str, Args ... args) { log(level::level_enum::debug, str, args...); } + template + void info( const char* str, Args ... args) { log(level::level_enum::info, str, args...); } + template + void warn( const char* str, Args ... args) { log(level::level_enum::warn, str, args...); } + template + void error( const char* str, Args ... args) { log(level::level_enum::err, str, args...); } + template + void critical( const char* str, Args ... args) { log(level::level_enum::critical, str, args...); } + }; + + + + + NGCORE_API std::shared_ptr GetLogger(const std::string& name); + NGCORE_API void SetLoggingLevel(level::level_enum level, const std::string& name); + NGCORE_API void AddFileSink(const std::string& filename, level::level_enum level, const std::string& logger); + NGCORE_API void AddConsoleSink(level::level_enum level, const std::string& logger); + NGCORE_API void ClearLoggingSinks(const std::string& logger); + NGCORE_API void FlushOnLoggingLevel(level::level_enum level, const std::string& logger); +} // namespace ngcore + +#endif // NETGEN_CORE_LOGGING_HPP diff --git a/libsrc/core/ngcore.hpp b/libsrc/core/ngcore.hpp index 1c8fa591..d73a3ed9 100644 --- a/libsrc/core/ngcore.hpp +++ b/libsrc/core/ngcore.hpp @@ -2,6 +2,10 @@ #define NETGEN_CORE_NGCORE_HPP #include "archive.hpp" +#include "exception.hpp" +#include "logging.hpp" +#include "profiler.hpp" +#include "symboltable.hpp" #include "version.hpp" #endif // NETGEN_CORE_NGCORE_HPP diff --git a/libsrc/core/ngcore_api.hpp b/libsrc/core/ngcore_api.hpp index c08d0ac9..a84a1381 100644 --- a/libsrc/core/ngcore_api.hpp +++ b/libsrc/core/ngcore_api.hpp @@ -5,8 +5,8 @@ #define NGCORE_API_EXPORT __declspec(dllexport) #define NGCORE_API_IMPORT __declspec(dllimport) #else - #define NGCORE_API_EXPORT - #define NGCORE_API_IMPORT + #define NGCORE_API_EXPORT __attribute__((visibility("default"))) + #define NGCORE_API_IMPORT __attribute__((visibility("default"))) #endif #ifdef NGCORE_EXPORTS @@ -15,15 +15,23 @@ #define NGCORE_API NGCORE_API_IMPORT #endif -namespace ngcore -{ -#if defined(__GNUC__) - inline bool likely (bool x) { return bool(__builtin_expect(long(x), 1L)); } - inline bool unlikely (bool x) { return bool(__builtin_expect(long(x), 0L)); } +#ifdef __INTEL_COMPILER + #ifdef WIN32 + #define NETGEN_INLINE __forceinline inline + #define NETGEN_LAMBDA_INLINE + #else + #define NETGEN_INLINE __forceinline inline + #define NETGEN_LAMBDA_INLINE __attribute__ ((__always_inline__)) + #endif #else - inline bool likely (bool x) { return x; } - inline bool unlikely (bool x) { return x; } + #ifdef __GNUC__ + #define NETGEN_INLINE __attribute__ ((__always_inline__)) inline + #define NETGEN_LAMBDA_INLINE __attribute__ ((__always_inline__)) + #define NETGEN_VLA + #else + #define NETGEN_INLINE inline + #define NETGEN_LAMBDA_INLINE + #endif #endif -} // namespace ngcore #endif // NETGEN_CORE_NGCORE_API_HPP diff --git a/libsrc/core/paje_trace.cpp b/libsrc/core/paje_trace.cpp new file mode 100644 index 00000000..844710d5 --- /dev/null +++ b/libsrc/core/paje_trace.cpp @@ -0,0 +1,667 @@ +#include +#include +#include +#include +#include +#include + +#include "archive.hpp" // for Demangle +#include "paje_trace.hpp" +#include "profiler.hpp" + +extern const char *header; + +namespace ngcore +{ + // Produce no traces by default + size_t PajeTrace::max_tracefile_size = 0; + + // If true, produce variable counting active threads + // increases trace by a factor of two + bool PajeTrace::trace_thread_counter = true; + bool PajeTrace::trace_threads = true; + + PajeTrace :: PajeTrace(int anthreads, std::string aname) + { + start_time = GetTimeCounter(); + + nthreads = anthreads; + tracefile_name = std::move(aname); + + int bytes_per_event=33; + max_num_events_per_thread = std::min( static_cast(std::numeric_limits::max()), max_tracefile_size/bytes_per_event/(2*nthreads+1)*10/7); + if(max_num_events_per_thread>0) + { + logger->info( "Tracefile size = {}MB", max_tracefile_size/1024/1024); + logger->info( "Tracing {} events per thread", max_num_events_per_thread); + } + + tasks.resize(nthreads); + int reserve_size = std::min(1000000U, max_num_events_per_thread); + for(auto & t : tasks) + t.reserve(reserve_size); + + links.resize(nthreads); + for(auto & l : links) + l.reserve(reserve_size); + + jobs.reserve(reserve_size); + timer_events.reserve(reserve_size); + + tracing_enabled = true; + } + + PajeTrace :: ~PajeTrace() + { + if(!tracefile_name.empty()) + Write(tracefile_name); + } + + + void PajeTrace::StopTracing() + { + if(tracing_enabled && max_num_events_per_thread>0) + { + logger->warn("Maximum number of traces reached, tracing is stopped now."); + } + tracing_enabled = false; + } + + class PajeFile + { + public: + static void Hue2RGB ( double x, double &r, double &g, double &b ) + { + double d = 1.0/6.0; + if(x logger = GetLogger("PajeTrace"); + + + double ConvertTime(TTimePoint t) { + // return time in milliseconds as double + // return std::chrono::duration(t-start_time).count()*1000.0; + // return std::chrono::duration(t-start_time).count() / 2.7e3; + return (t-start_time) / 2.7e6; + } + + enum PType + { + SET_VARIABLE=1, + ADD_VARIABLE, + SUB_VARIABLE, + PUSH_STATE, + POP_STATE, + START_LINK, + STOP_LINK + }; + + struct PajeEvent + { + PajeEvent( int aevent_type, double atime, int atype, int acontainer, double avar_value ) + : time(atime), var_value(avar_value), event_type(aevent_type), type(atype), container(acontainer) + { } + + PajeEvent( int aevent_type, double atime, int atype, int acontainer, int avalue = 0, int aid = 0, bool avalue_is_alias = true ) + : time(atime), event_type(aevent_type), type(atype), container(acontainer), value(avalue), id(aid), value_is_alias(avalue_is_alias) + { } + + PajeEvent( int aevent_type, double atime, int atype, int acontainer, int avalue, int astart_container, int akey ) + : time(atime), event_type(aevent_type), type(atype), container(acontainer), value(avalue), start_container(astart_container), id(akey) + { } + + double time; + double var_value = 0.0; + int event_type; + int type; + int container; + int value = 0; + int start_container = 0; + int id = 0; + bool value_is_alias = true; + + bool operator < (const PajeEvent & other) const { + // Same start and stop times can occur for very small tasks -> take "starting" events first (eg. PajePushState before PajePopState) + if(time == other.time) + return event_type < other.event_type; + return (time < other.time); + } + + int write(FILE *stream) + { + const int &key = id; + const int &end_container = start_container; + switch(event_type) + { + case PajeSetVariable: + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\t%.15g\n", PajeSetVariable, time, type, container, var_value ); // NOLINT + case PajeAddVariable: + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\t%.15g\n", PajeAddVariable, time, type, container, var_value ); // NOLINT + case PajeSubVariable: + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\t%.15g\n", PajeSubVariable, time, type, container, var_value ); // NOLINT + case PajePushState: + if(value_is_alias) + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\ta%d\t%d\n", PajePushState, time, type, container, value, id); // NOLINT + else + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\t%d\t%d\n", PajePushState, time, type, container, value, id); // NOLINT + case PajePopState: + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\n", PajePopState, time, type, container ); // NOLINT + case PajeStartLink: + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\t%d\ta%d\t%d\n", PajeStartLink, time, type, container, value, start_container, key ); // NOLINT + case PajeEndLink: + return fprintf( stream, "%d\t%.15g\ta%d\ta%d\t%d\ta%d\t%d\n", PajeEndLink, time, type, container, value, end_container, key ); // NOLINT + } + return 0; + } + }; + + std::vector events; + + public: + PajeFile() = delete; + PajeFile(const PajeFile &) = delete; + PajeFile(PajeFile &&) = delete; + void operator=(const PajeFile &) = delete; + void operator=(PajeFile &&) = delete; + + PajeFile( const std::string & filename, TTimePoint astart_time ) + { + start_time = astart_time; + ctrace_stream = fopen (filename.c_str(),"w"); // NOLINT + fprintf(ctrace_stream, "%s", header ); // NOLINT + alias_counter = 0; + } + + ~PajeFile() + { + fclose (ctrace_stream); // NOLINT + } + + int DefineContainerType ( int parent_type, const std::string & name ) + { + int alias = ++alias_counter; + if(parent_type!=0) + fprintf( ctrace_stream, "%d\ta%d\ta%d\t\"%s\"\n", PajeDefineContainerType, alias, parent_type, name.c_str() ); // NOLINT + else + fprintf( ctrace_stream, "%d\ta%d\t%d\t\"%s\"\n", PajeDefineContainerType, alias, parent_type, name.c_str() ); // NOLINT + return alias; + } + + int DefineVariableType ( int container_type, const std::string & name ) + { + int alias = ++alias_counter; + fprintf( ctrace_stream, "%d\ta%d\ta%d\t\"%s\"\t\"1.0 1.0 1.0\"\n", PajeDefineVariableType, alias, container_type, name.c_str() ); // NOLINT + return alias; + } + + int DefineStateType ( int type, const std::string & name ) + { + int alias = ++alias_counter; + fprintf( ctrace_stream, "%d\ta%d\ta%d\t\"%s\"\n", PajeDefineStateType, alias, type, name.c_str() ); // NOLINT + return alias; + } + + // int DefineEventType () + // { + // Write("event not implemented"); + // } + + int DefineLinkType (int parent_container_type, int start_container_type, int stop_container_type, const std::string & name) + { + int alias = ++alias_counter; + fprintf( ctrace_stream, "%d\ta%d\ta%d\ta%d\ta%d\t\"%s\"\n", PajeDefineLinkType, alias, parent_container_type, start_container_type, stop_container_type, name.c_str() ); // NOLINT + return alias; + } + + int DefineEntityValue (int type, const std::string & name, double hue = -1) + { + if(hue==-1) + { + std::hash shash; + size_t h = shash(name); + h ^= h>>32U; + h = static_cast(h); + hue = h*1.0/std::numeric_limits::max(); + } + + int alias = ++alias_counter; + double r,g,b; + Hue2RGB( hue, r, g, b ); + fprintf( ctrace_stream, "%d\ta%d\ta%d\t\"%s\"\t\"%.15g %.15g %.15g\"\n", PajeDefineEntityValue, alias, type, name.c_str(), r,g,b ); // NOLINT + return alias; + } + + int CreateContainer ( int type, int parent, const std::string & name ) + { + int alias = ++alias_counter; + if(parent!=0) + fprintf( ctrace_stream, "%d\t0\ta%d\ta%d\ta%d\t\"%s\"\n", PajeCreateContainer, alias, type, parent, name.c_str() ); // NOLINT + else + fprintf( ctrace_stream, "%d\t0\ta%d\ta%d\t%d\t\"%s\"\n", PajeCreateContainer, alias, type, parent, name.c_str() ); // NOLINT + return alias; + } + void DestroyContainer () + {} + + void SetVariable (TTimePoint time, int type, int container, double value ) + { + events.emplace_back( PajeEvent( PajeSetVariable, ConvertTime(time), type, container, value ) ); + } + + void AddVariable (TTimePoint time, int type, int container, double value ) + { + events.emplace_back( PajeEvent( PajeAddVariable, ConvertTime(time), type, container, value ) ); + } + + void SubVariable (TTimePoint time, int type, int container, double value ) + { + events.emplace_back( PajeEvent( PajeSubVariable, ConvertTime(time), type, container, value ) ); + } + + void SetState () + {} + + void PushState ( TTimePoint time, int type, int container, int value, int id = 0, bool value_is_alias = true ) + { + events.emplace_back( PajeEvent( PajePushState, ConvertTime(time), type, container, value, id, value_is_alias) ); + } + + void PopState ( TTimePoint time, int type, int container ) + { + events.emplace_back( PajeEvent( PajePopState, ConvertTime(time), type, container ) ); + } + + void ResetState () + {} + + void StartLink ( TTimePoint time, int type, int container, int value, int start_container, int key ) + { + events.emplace_back( PajeEvent( PajeStartLink, ConvertTime(time), type, container, value, start_container, key ) ); + } + + void EndLink ( TTimePoint time, int type, int container, int value, int end_container, int key ) + { + events.emplace_back( PajeEvent( PajeEndLink, ConvertTime(time), type, container, value, end_container, key ) ); + } + + void NewEvent () + {} + + void WriteEvents() + { + logger->info("Sorting traces..."); + std::sort (events.begin(), events.end()); + + logger->info("Writing traces... "); + for (auto & event : events) + { + event.write( ctrace_stream ); +// fprintf( ctrace_stream, "%s", buf ); // NOLINT + } + logger->info("Done"); + } + + private: + enum + { + PajeDefineContainerType = 0, + PajeDefineVariableType = 1, + PajeDefineStateType = 2, + PajeDefineEventType = 3, + PajeDefineLinkType = 4, + PajeDefineEntityValue = 5, + PajeCreateContainer = 6, + PajeDestroyContainer = 7, + PajeSetVariable = 8, + PajeAddVariable = 9, + PajeSubVariable = 10, + PajeSetState = 11, + PajePushState = 12, + PajePopState = 13, + PajeResetState = 14, + PajeStartLink = 15, + PajeEndLink = 16, + PajeNewEvent = 17 + }; + + }; + + NGCORE_API PajeTrace *trace; + + void PajeTrace::Write( const std::string & filename ) + { + int n_events = jobs.size() + timer_events.size(); + for(auto & vtasks : tasks) + n_events += vtasks.size(); + + logger->info("{} events traced", n_events); + + if(n_events==0) + { + logger->info("No data traced, skip writing trace file"); + return; + } + + if(!tracing_enabled) + { + logger->warn("Tracing stopped during computation due to tracefile size limit of {} megabytes.", max_tracefile_size/1024/1024); + } + + PajeFile paje(filename, start_time); + + const int container_type_task_manager = paje.DefineContainerType( 0, "Task Manager" ); + const int container_type_node = paje.DefineContainerType( container_type_task_manager, "Node"); + const int container_type_thread = paje.DefineContainerType( container_type_task_manager, "Thread"); + const int container_type_timer = container_type_thread; //paje.DefineContainerType( container_type_task_manager, "Timers"); + const int container_type_jobs = paje.DefineContainerType( container_type_task_manager, "Jobs"); + + const int state_type_job = paje.DefineStateType( container_type_jobs, "Job" ); + const int state_type_task = paje.DefineStateType( container_type_thread, "Task" ); + const int state_type_timer = paje.DefineStateType( container_type_timer, "Timer state" ); + + const int variable_type_active_threads = paje.DefineVariableType( container_type_jobs, "Active threads" ); + + const int container_task_manager = paje.CreateContainer( container_type_task_manager, 0, "The task manager" ); + const int container_jobs = paje.CreateContainer( container_type_jobs, container_task_manager, "Jobs" ); + paje.SetVariable( start_time, variable_type_active_threads, container_jobs, 0.0 ); + + const int num_nodes = 1; //task_manager ? task_manager->GetNumNodes() : 1; + + std::vector container_nodes; + container_nodes.reserve(num_nodes); + for(int i=0; i thread_aliases; + thread_aliases.reserve(nthreads); + if(trace_threads) + for (int i=0; i job_map; + std::map job_task_map; + + for(Job & j : jobs) + if(job_map.find(j.type) == job_map.end()) + { + std::string name = Demangle(j.type->name()); + job_map[j.type] = paje.DefineEntityValue( state_type_job, name, -1 ); + job_task_map[j.type] = paje.DefineEntityValue( state_type_task, name, -1 ); + } + + for(Job & j : jobs) + { + paje.PushState( j.start_time, state_type_job, container_jobs, job_map[j.type] ); + paje.PopState( j.stop_time, state_type_job, container_jobs ); + } + + std::set timer_ids; + std::map timer_aliases; + + for(auto & event : timer_events) + timer_ids.insert(event.timer_id); + + + for(auto & vtasks : tasks) + for (Task & t : vtasks) + if(t.id_type==Task::ID_TIMER) + timer_ids.insert(t.id); + + for(auto id : timer_ids) + timer_aliases[id] = paje.DefineEntityValue( state_type_timer, NgProfiler::GetName(id), -1 ); + + int timerdepth = 0; + int maxdepth = 0; + for(auto & event : timer_events) + { + if(event.is_start) + { + timerdepth++; + maxdepth = timerdepth>maxdepth ? timerdepth : maxdepth; + } + else + timerdepth--; + } + + std::vector timer_container_aliases; + timer_container_aliases.resize(maxdepth); + for(int i=0; i links_merged; + links_merged.reserve(nlinks); + std::vector pos(nthreads); + + int nlinks_merged = 0; + while(nlinks_merged < nlinks) + { + int minpos = -1; + TTimePoint mintime = -1; + for (int t = 0; t started_links; + + int link_type = paje.DefineLinkType(container_type_node, container_type_thread, container_type_thread, "links"); + + // match links + for ( auto & l : links_merged ) + { + if(l.is_start) + { + started_links.push_back(l); + } + else + { + unsigned int i = 0; + while(i +#include + +#include "logging.hpp" // for logger +#include "ngcore_api.hpp" // for NGCORE_API +#include "utils.hpp" + +namespace ngcore +{ + + extern NGCORE_API class PajeTrace *trace; + class PajeTrace + { + public: + using TClock = std::chrono::system_clock; + + protected: + std::shared_ptr logger = GetLogger("PajeTrace"); + private: + NGCORE_API static size_t max_tracefile_size; + NGCORE_API static bool trace_thread_counter; + NGCORE_API static bool trace_threads; + + bool tracing_enabled; + TTimePoint start_time; + int nthreads; + + public: + + // Approximate number of events to trace. Tracing will + // be stopped if any thread reaches this number of events + unsigned int max_num_events_per_thread; + + static void SetTraceThreads( bool atrace_threads ) + { + trace_threads = atrace_threads; + } + + static void SetTraceThreadCounter( bool trace_threads ) + { + trace_thread_counter = trace_threads; + } + + static void SetMaxTracefileSize( size_t max_size ) + { + max_tracefile_size = max_size; + } + + std::string tracefile_name; + + struct Job + { + int job_id; + const std::type_info *type; + TTimePoint start_time; + TTimePoint stop_time; + }; + + struct Task + { + int thread_id; + + int id; + int id_type; + + int additional_value; + + TTimePoint start_time; + TTimePoint stop_time; + + static constexpr int ID_NONE = -1; + static constexpr int ID_JOB = 1; + static constexpr int ID_TIMER = 2; + }; + + struct TimerEvent + { + int timer_id; + TTimePoint time; + bool is_start; + int thread_id; + + bool operator < (const TimerEvent & other) const { return time < other.time; } + }; + + struct ThreadLink + { + int thread_id; + int key; + TTimePoint time; + bool is_start; + bool operator < (const ThreadLink & other) const { return time < other.time; } + }; + + std::vector > tasks; + std::vector jobs; + std::vector timer_events; + std::vector > links; + + public: + NGCORE_API void StopTracing(); + + PajeTrace() = delete; + PajeTrace(const PajeTrace &) = delete; + PajeTrace(PajeTrace &&) = delete; + NGCORE_API PajeTrace(int anthreads, std::string aname = ""); + NGCORE_API ~PajeTrace(); + + void operator=(const PajeTrace &) = delete; + void operator=(PajeTrace &&) = delete; + + void StartTimer(int timer_id) + { + if(!tracing_enabled) return; + if(unlikely(timer_events.size() == max_num_events_per_thread)) + StopTracing(); + timer_events.push_back(TimerEvent{timer_id, GetTimeCounter(), true}); + } + + void StopTimer(int timer_id) + { + if(!tracing_enabled) return; + if(unlikely(timer_events.size() == max_num_events_per_thread)) + StopTracing(); + timer_events.push_back(TimerEvent{timer_id, GetTimeCounter(), false}); + } + + NETGEN_INLINE int StartTask(int thread_id, int id, int id_type = Task::ID_NONE, int additional_value = -1) + { + if(!tracing_enabled) return -1; + if(!trace_threads && !trace_thread_counter) return -1; + if(unlikely(tasks[thread_id].size() == max_num_events_per_thread)) + StopTracing(); + int task_num = tasks[thread_id].size(); + tasks[thread_id].push_back( Task{thread_id, id, id_type, additional_value, GetTimeCounter()} ); + return task_num; + } + + void StopTask(int thread_id, int task_num) + { + if(!trace_threads && !trace_thread_counter) return; + if(task_num>=0) + tasks[thread_id][task_num].stop_time = GetTimeCounter(); + } + + void SetTask(int thread_id, int task_num, int additional_value) { + if(!trace_threads && !trace_thread_counter) return; + if(task_num>=0) + tasks[thread_id][task_num].additional_value = additional_value; + } + + void StartJob(int job_id, const std::type_info & type) + { + if(!tracing_enabled) return; + if(jobs.size() == max_num_events_per_thread) + StopTracing(); + jobs.push_back( Job{job_id, &type, GetTimeCounter()} ); + } + + void StopJob() + { + if(tracing_enabled) + jobs.back().stop_time = GetTimeCounter(); + } + + void StartLink(int thread_id, int key) + { + if(!tracing_enabled) return; + if(links[thread_id].size() == max_num_events_per_thread) + StopTracing(); + links[thread_id].push_back( ThreadLink{thread_id, key, GetTimeCounter(), true} ); + } + + void StopLink(int thread_id, int key) + { + if(!tracing_enabled) return; + if(links[thread_id].size() == max_num_events_per_thread) + StopTracing(); + links[thread_id].push_back( ThreadLink{thread_id, key, GetTimeCounter(), false} ); + } + + void Write( const std::string & filename ); + + }; +} // namespace ngcore + +#endif // NETGEN_CORE_PAJE_TRACE_HPP diff --git a/libsrc/core/profiler.cpp b/libsrc/core/profiler.cpp new file mode 100644 index 00000000..3302342d --- /dev/null +++ b/libsrc/core/profiler.cpp @@ -0,0 +1,117 @@ +#include + +#include "profiler.hpp" + +namespace ngcore +{ + std::vector NgProfiler::timers(NgProfiler::SIZE); // NOLINT + + std::string NgProfiler::filename; + + size_t NgProfiler::dummy_thread_times[NgProfiler::SIZE]; + size_t * NgProfiler::thread_times = NgProfiler::dummy_thread_times; // NOLINT + size_t NgProfiler::dummy_thread_flops[NgProfiler::SIZE]; + size_t * NgProfiler::thread_flops = NgProfiler::dummy_thread_flops; // NOLINT + + std::shared_ptr NgProfiler::logger = GetLogger("Profiler"); // NOLINT + + NgProfiler :: NgProfiler() + { + for (auto & t : timers) + { + t.tottime = 0.0; + t.usedcounter = 0; + t.flops = 0.0; + } + } + + NgProfiler :: ~NgProfiler() + { + if (filename.length()) + { + logger->debug( "write profile to file {}", filename ); + FILE *prof = fopen(filename.c_str(),"w"); // NOLINT + Print (prof); + fclose(prof); // NOLINT + } + + if (getenv ("NGPROFILE")) + { + std::string filename = "netgen.prof"; +#ifdef PARALLEL + filename += "."+ToString(id); +#endif + if (id == 0) logger->info( "write profile to file {}", filename ); + FILE *prof = fopen(filename.c_str(),"w"); // NOLINT + Print (prof); + fclose(prof); // NOLINT + } + } + + void NgProfiler :: Print (FILE * prof) + { + int i = 0; + for (auto & t : timers) + { + if (t.count != 0 || t.usedcounter != 0) + { + fprintf(prof,"job %3i calls %8li, time %6.4f sec",i,t.count,t.tottime); // NOLINT + if(t.flops != 0.0) + fprintf(prof,", MFlops = %6.2f",t.flops / (t.tottime) * 1e-6); // NOLINT + if(t.loads != 0.0) + fprintf(prof,", MLoads = %6.2f",t.loads / (t.tottime) * 1e-6); // NOLINT + if(t.stores != 0.0) + fprintf(prof,", MStores = %6.2f",t.stores / (t.tottime) * 1e-6); // NOLINT + if(t.usedcounter) + fprintf(prof," %s",t.name.c_str()); // NOLINT + fprintf(prof,"\n"); // NOLINT + } + i++; + } + } + + + int NgProfiler :: CreateTimer (const std::string & name) + { + static std::mutex createtimer_mutex; + int nr = -1; + { + std::lock_guard guard(createtimer_mutex); + for (int i = SIZE-1; i > 0; i--) + { + auto & t = timers[i]; + if (!t.usedcounter) + { + t.usedcounter = 1; + t.name = name; + nr = i; + break; + } + } + } + if (nr > -1) return nr; + static bool first_overflow = true; + if (first_overflow) + { + first_overflow = false; + NgProfiler::logger->warn("no more timer available, reusing last one"); + } + return 0; + } + + void NgProfiler :: Reset () + { + for(auto & t : timers) + { + t.tottime = 0.0; + t.count = 0; + t.flops = 0.0; + t.loads = 0; + t.stores = 0; + } + } + + NgProfiler prof; // NOLINT + + +} // namespace ngcore diff --git a/libsrc/core/profiler.hpp b/libsrc/core/profiler.hpp new file mode 100644 index 00000000..5a9b16e1 --- /dev/null +++ b/libsrc/core/profiler.hpp @@ -0,0 +1,304 @@ +#ifndef NETGEN_CORE_PROFILER_HPP +#define NETGEN_CORE_PROFILER_HPP + +#include +#include + +#include "logging.hpp" +#include "paje_trace.hpp" +#include "utils.hpp" + +namespace ngcore +{ + class NgProfiler + { + public: + /// maximal number of timers + enum { SIZE = 8*1024 }; + + struct TimerVal + { + TimerVal() = default; + + double tottime = 0.0; + double starttime = 0.0; + double flops = 0.0; + double loads = 0.0; + double stores = 0.0; + long count = 0; + std::string name = ""; + int usedcounter = 0; + }; + + NGCORE_API static std::vector timers; + + NGCORE_API static TTimePoint * thread_times; + NGCORE_API static TTimePoint * thread_flops; + NGCORE_API static std::shared_ptr logger; + NGCORE_API static size_t dummy_thread_times[NgProfiler::SIZE]; + NGCORE_API static size_t dummy_thread_flops[NgProfiler::SIZE]; + private: + + NGCORE_API static std::string filename; + public: + NgProfiler(); + ~NgProfiler(); + + NgProfiler(const NgProfiler &) = delete; + NgProfiler(NgProfiler &&) = delete; + void operator=(const NgProfiler &) = delete; + void operator=(NgProfiler &&) = delete; + + static void SetFileName (const std::string & afilename) { filename = afilename; } + + /// create new timer, use integer index + NGCORE_API static int CreateTimer (const std::string & name); + + NGCORE_API static void Reset (); + + + /// start timer of index nr + static void StartTimer (int nr) + { + timers[nr].starttime = WallTime(); timers[nr].count++; + } + + /// stop timer of index nr + static void StopTimer (int nr) + { + timers[nr].tottime += WallTime()-timers[nr].starttime; + } + + static void StartThreadTimer (size_t nr, size_t tid) + { + thread_times[tid*SIZE+nr] -= GetTimeCounter(); // NOLINT + } + + static void StopThreadTimer (size_t nr, size_t tid) + { + thread_times[tid*SIZE+nr] += GetTimeCounter(); // NOLINT + } + + static void AddThreadFlops (size_t nr, size_t tid, size_t flops) + { + thread_flops[tid*SIZE+nr] += flops; // NOLINT + } + + /// if you know number of flops, provide them to obtain the MFlop - rate + static void AddFlops (int nr, double aflops) { timers[nr].flops += aflops; } + static void AddLoads (int nr, double aloads) { timers[nr].loads += aloads; } + static void AddStores (int nr, double astores) { timers[nr].stores += astores; } + + static int GetNr (const std::string & name) + { + for (int i = SIZE-1; i >= 0; i--) + if (timers[i].name == name) + return i; + return -1; + } + + static double GetTime (int nr) + { + return timers[nr].tottime; + } + + static double GetTime (const std::string & name) + { + for (int i = SIZE-1; i >= 0; i--) + if (timers[i].name == name) + return GetTime (i); + return 0; + } + + static long int GetCounts (int nr) + { + return timers[nr].count; + } + + static double GetFlops (int nr) + { + return timers[nr].flops; + } + + /// change name + static void SetName (int nr, const std::string & name) { timers[nr].name = name; } + static std::string GetName (int nr) { return timers[nr].name; } + /// print profile + NGCORE_API static void Print (FILE * prof); + + class RegionTimer + { + int nr; + public: + /// start timer + RegionTimer (int anr) : nr(anr) { NgProfiler::StartTimer(nr); } + /// stop timer + ~RegionTimer () { NgProfiler::StopTimer(nr); } + + RegionTimer() = delete; + RegionTimer(const RegionTimer &) = delete; + RegionTimer(RegionTimer &&) = delete; + void operator=(const RegionTimer &) = delete; + void operator=(RegionTimer &&) = delete; + }; + }; + + + + class NGCORE_API Timer + { + int timernr; + int priority; + public: + Timer (const std::string & name, int apriority = 1) + : priority(apriority) + { + timernr = NgProfiler::CreateTimer (name); + } + void SetName (const std::string & name) + { + NgProfiler::SetName (timernr, name); + } + void Start () + { + if (priority <= 2) + NgProfiler::StartTimer (timernr); + if (priority <= 1) + if(trace) trace->StartTimer(timernr); + } + void Stop () + { + if (priority <= 2) + NgProfiler::StopTimer (timernr); + if (priority <= 1) + if(trace) trace->StopTimer(timernr); + } + void AddFlops (double aflops) + { + if (priority <= 2) + NgProfiler::AddFlops (timernr, aflops); + } + + double GetTime () { return NgProfiler::GetTime(timernr); } + long int GetCounts () { return NgProfiler::GetCounts(timernr); } + double GetMFlops () + { return NgProfiler::GetFlops(timernr) + / NgProfiler::GetTime(timernr) * 1e-6; } + operator int () { return timernr; } + }; + + + /** + Timer object. + Start / stop timer at constructor / destructor. + */ + class RegionTimer + { + Timer & timer; + public: + /// start timer + RegionTimer (Timer & atimer) : timer(atimer) { timer.Start(); } + /// stop timer + ~RegionTimer () { timer.Stop(); } + + RegionTimer() = delete; + RegionTimer(const RegionTimer &) = delete; + RegionTimer(RegionTimer &&) = delete; + void operator=(const RegionTimer &) = delete; + void operator=(RegionTimer &&) = delete; + }; + + class ThreadRegionTimer + { + size_t nr; + size_t tid; + public: + /// start timer + ThreadRegionTimer (size_t _nr, size_t _tid) : nr(_nr), tid(_tid) + { NgProfiler::StartThreadTimer(nr, tid); } + /// stop timer + ~ThreadRegionTimer () + { NgProfiler::StopThreadTimer(nr, tid); } + + ThreadRegionTimer() = delete; + ThreadRegionTimer(ThreadRegionTimer &&) = delete; + ThreadRegionTimer(const ThreadRegionTimer &) = delete; + void operator=(const ThreadRegionTimer &) = delete; + void operator=(ThreadRegionTimer &&) = delete; + }; + + class RegionTracer + { + int nr; + int thread_id; + public: + static constexpr int ID_JOB = PajeTrace::Task::ID_JOB; + static constexpr int ID_NONE = PajeTrace::Task::ID_NONE; + static constexpr int ID_TIMER = PajeTrace::Task::ID_TIMER; + + RegionTracer() = delete; + RegionTracer(RegionTracer &&) = delete; + RegionTracer(const RegionTracer &) = delete; + void operator=(const RegionTracer &) = delete; + void operator=(RegionTracer &&) = delete; + + /// start trace + RegionTracer (int athread_id, int region_id, int id_type = ID_NONE, int additional_value = -1 ) + : thread_id(athread_id) + { + if (trace) + nr = trace->StartTask (athread_id, region_id, id_type, additional_value); + } + /// start trace with timer + RegionTracer (int athread_id, Timer & timer, int additional_value = -1 ) + : thread_id(athread_id) + { + if (trace) + nr = trace->StartTask (athread_id, static_cast(timer), ID_TIMER, additional_value); + } + + /// set user defined value + void SetValue( int additional_value ) + { + if (trace) + trace->SetTask( thread_id, nr, additional_value ); + } + + /// stop trace + ~RegionTracer () + { + if (trace) + trace->StopTask (thread_id, nr); + } + }; + + + // Helper function for timings + // Run f() at least min_iterations times until max_time seconds elapsed + // returns minimum runtime for a call of f() + template + double RunTiming( TFunc f, double max_time = 0.5, int min_iterations = 10 ) + { + // Make sure the whole test run does not exceed maxtime + double tend = WallTime()+max_time; + + // warmup + f(); + + double tres = std::numeric_limits::max(); + int iteration = 0; + while(WallTime() + +#include "logging.hpp" + +namespace py = pybind11; +using namespace ngcore; + +PYBIND11_MODULE(pyngcore, m) // NOLINT +{ + py::enum_(m, "LOG_LEVEL", "Logging level") + .value("Trace", level::trace) + .value("Debug", level::debug) + .value("Info", level::info) + .value("Warn", level::warn) + .value("Error", level::err) + .value("Critical", level::critical) + .value("Off", level::off); + + m.def("SetLoggingLevel", &SetLoggingLevel, py::arg("level"), py::arg("logger")="", + "Set logging level, if name is given only to the specific logger, else set the global logging level"); + m.def("AddFileSink", &AddFileSink, py::arg("filename"), py::arg("level"), py::arg("logger")="", + "Add File sink, either only to logger specified or globally to all loggers"); + m.def("AddConsoleSink", &AddConsoleSink, py::arg("level"), py::arg("logger")="", + "Add console output for specific logger or all if none given"); + m.def("ClearLoggingSinks", &ClearLoggingSinks, py::arg("logger")="", + "Clear sinks of specific logger, or all if none given"); + m.def("FlushOnLoggingLevel", &FlushOnLoggingLevel, py::arg("level"), py::arg("logger")="", + "Flush every message with level at least `level` for specific logger or all loggers if none given."); +} diff --git a/libsrc/core/symboltable.hpp b/libsrc/core/symboltable.hpp new file mode 100644 index 00000000..fdab713c --- /dev/null +++ b/libsrc/core/symboltable.hpp @@ -0,0 +1,144 @@ +#ifndef NETGEN_CORE_SYMBOLTABLE_HPP +#define NETGEN_CORE_SYMBOLTABLE_HPP + +#include +#include +#include + +#include "archive.hpp" +#include "exception.hpp" +#include "ngcore_api.hpp" + +namespace ngcore +{ + /** + A symbol table. + + The symboltable provides a mapping from string identifiers + to the generic type T. The strings are copied. + Complexity by name access is linear, by index is constant. + */ + template + class SymbolTable + { + std::vector names; + std::vector data; + public: + using value_type = T; + using reference = typename std::vector::reference; + using const_reference = typename std::vector::const_reference; + + /// Creates a symboltable + SymbolTable () = default; + SymbolTable (const SymbolTable &) = default; + SymbolTable (SymbolTable &&) noexcept = default; + + ~SymbolTable() = default; + + SymbolTable& operator=(const SymbolTable&) = default; + SymbolTable& operator=(SymbolTable&&) = default; + + template + auto DoArchive(Archive& ar) -> typename std::enable_if, void>::type + { + ar & names & data; + } + + /// INDEX of symbol name, throws exception if unused + size_t Index (const std::string & name) const + { + for (size_t i = 0; i < names.size(); i++) + if (names[i] == name) return i; + throw RangeException("SymbolTable", name); + } + + /// Index of symbol name, returns -1 if unused + int CheckIndex (const std::string & name) const + { + for (int i = 0; i < names.size(); i++) + if (names[i] == name) return i; + return -1; + } + + /// number of identifiers + size_t Size() const + { + return data.size(); + } + + /// Returns reference to element. exception for unused identifier + reference operator[] (const std::string & name) + { + return data[Index (name)]; + } + + const_reference operator[] (const std::string & name) const + { + return data[Index (name)]; + } + + /// Returns reference to i-th element, range check only in debug build + reference operator[] (size_t i) + { + NETGEN_CHECK_RANGE(i, 0, data.size()); + return data[i]; + } + + /// Returns const reference to i-th element, range check only in debug build + const_reference operator[] (size_t i) const + { + NETGEN_CHECK_RANGE(i, 0, data.size()); + return data[i]; + } + + /// Returns name of i-th element, range check only in debug build + const std::string & GetName (size_t i) const + { + NETGEN_CHECK_RANGE(i, 0, names.size()); + return names[i]; + } + + /// Associates el to the string name, overrides if name is used + void Set (const std::string & name, const T & el) + { + int i = CheckIndex (name); + if (i >= 0) + data[i] = el; + else + { + data.push_back(el); + names.push_back(name); + } + } + + bool Used (const std::string & name) const + { + return CheckIndex(name) >= 0; + } + + /// Deletes symboltable + inline void DeleteAll () + { + names.clear(); + data.clear(); + } + + // Adds all elements from other symboltable + SymbolTable& Update(const SymbolTable& tbl2) + { + for (size_t i = 0; i < tbl2.Size(); i++) + Set (tbl2.GetName(i), tbl2[i]); + return *this; + } + }; + + template + std::ostream & operator<< (std::ostream & ost, const SymbolTable & st) + { + for (int i = 0; i < st.Size(); i++) + ost << st.GetName(i) << " : " << st[i] << std::endl; + return ost; + } +} // namespace ngcore + +#endif // NETGEN_CORE_SYMBOLTABLE_HPP diff --git a/libsrc/core/utils.cpp b/libsrc/core/utils.cpp new file mode 100644 index 00000000..0347d244 --- /dev/null +++ b/libsrc/core/utils.cpp @@ -0,0 +1,42 @@ +#include "utils.hpp" +#include "logging.hpp" + +#ifndef WIN32 +#include +#endif +#include + +namespace ngcore +{ + // parallel netgen + int id = 0, ntasks = 1; + +#ifdef WIN32 + // windows does demangling in typeid(T).name() + NGCORE_API std::string Demangle(const char* typeinfo) { return typeinfo; } +#else + NGCORE_API std::string Demangle(const char* typeinfo) { int status; return abi::__cxa_demangle(typeinfo, + nullptr, + nullptr, + &status); } +#endif + + double ticks_per_second = [] () noexcept + { + auto tick_start = GetTimeCounter(); + double tstart = WallTime(); + double tend = WallTime()+0.001; + + // wait for 1ms and compare wall time with time counter + while(WallTime() wall_time_start = TClock::now(); + +} // namespace ngcore + diff --git a/libsrc/core/utils.hpp b/libsrc/core/utils.hpp new file mode 100644 index 00000000..f2a73578 --- /dev/null +++ b/libsrc/core/utils.hpp @@ -0,0 +1,70 @@ +#ifndef NETGEN_CORE_UTILS_HPP +#define NETGEN_CORE_UTILS_HPP + +#include +#include +#include +#include +#include + +#ifdef WIN32 +#include // for __rdtsc() CPU time step counter +#else +#include // for __rdtsc() CPU time step counter +#endif // WIN32 + +#include "ngcore_api.hpp" // for NGCORE_API + +namespace ngcore +{ + // MPI rank, nranks TODO: Rename + extern NGCORE_API int id, ntasks; + + NGCORE_API std::string Demangle(const char* typeinfo); + +#if defined(__GNUC__) + inline bool likely (bool x) { return bool(__builtin_expect(long(x), 1L)); } + inline bool unlikely (bool x) { return bool(__builtin_expect(long(x), 0L)); } +#else + inline bool likely (bool x) { return x; } + inline bool unlikely (bool x) { return x; } +#endif + + using TClock = std::chrono::system_clock; + extern NGCORE_API const std::chrono::time_point wall_time_start; + + // Time in seconds since program start + inline double WallTime () noexcept + { + std::chrono::time_point now = TClock::now(); + std::chrono::duration elapsed_seconds = now-wall_time_start; + return elapsed_seconds.count(); + } + + // High precision clock counter register + using TTimePoint = size_t; + extern NGCORE_API double ticks_per_second; + + inline TTimePoint GetTimeCounter() noexcept + { + return TTimePoint(__rdtsc()); + } + + template + inline std::string ToString (const T& t) + { + std::stringstream ss; + ss << t; + return ss.str(); + } + + template + std::ostream& operator << (std::ostream& ost, const std::map& map) + { + for(auto& val : map) + ost << "\n" << val.first << ": " << val.second; + return ost; + } +} // namespace ngcore + +#endif // NETGEN_CORE_UTILS_HPP diff --git a/libsrc/core/version.hpp b/libsrc/core/version.hpp index 69598716..aea50bf6 100644 --- a/libsrc/core/version.hpp +++ b/libsrc/core/version.hpp @@ -1,6 +1,7 @@ #ifndef NETGEN_CORE_VERSION_HPP #define NETGEN_CORE_VERSION_HPP +#include #include #include @@ -83,6 +84,11 @@ namespace ngcore bool operator <=(const VersionInfo& other) const { return !((*this) > other); } bool operator >=(const VersionInfo& other) const { return !((*this) < other); } }; + + inline std::ostream& operator << (std::ostream& ost, const VersionInfo& version) + { + return ost << version.to_string(); + } } // namespace ngcore #endif // NETGEN_CORE_VERSION_HPP diff --git a/libsrc/csg/CMakeLists.txt b/libsrc/csg/CMakeLists.txt index 224a48e4..eaf7dffc 100644 --- a/libsrc/csg/CMakeLists.txt +++ b/libsrc/csg/CMakeLists.txt @@ -11,12 +11,10 @@ if(APPLE) set_target_properties( csg PROPERTIES SUFFIX ".so") endif(APPLE) -if(NOT WIN32) - target_link_libraries(csg mesh ${PYTHON_LIBRARIES}) - target_link_libraries(csg ${PYTHON_LIBRARIES}) - install( TARGETS csg ${NG_INSTALL_DIR}) -endif(NOT WIN32) +target_link_libraries(csg PUBLIC mesh ${PYTHON_LIBRARIES}) +install( TARGETS csg ${NG_INSTALL_DIR}) +target_link_libraries(csg PUBLIC ngcore) if(USE_GUI) add_library(csgvis ${NG_LIB_TYPE} vscsg.cpp ) diff --git a/libsrc/csg/csgeom.cpp b/libsrc/csg/csgeom.cpp index 87cae1b4..4aebe271 100644 --- a/libsrc/csg/csgeom.cpp +++ b/libsrc/csg/csgeom.cpp @@ -562,7 +562,7 @@ namespace netgen const Surface * CSGeometry :: GetSurface (const char * name) const { if (surfaces.Used(name)) - return surfaces.Get(name); + return surfaces[name]; else return NULL; } @@ -585,7 +585,7 @@ namespace netgen Solid * oldsol = NULL; if (solids.Used (name)) - oldsol = solids.Get(name); + oldsol = solids[name]; solids.Set (name, sol); sol->SetName (name); @@ -605,7 +605,7 @@ namespace netgen const Solid * CSGeometry :: GetSolid (const char * name) const { if (solids.Used(name)) - return solids.Get(name); + return solids[name]; else return NULL; } @@ -616,8 +616,8 @@ namespace netgen const Solid * CSGeometry :: GetSolid (const string & name) const { - if (solids.Used(name.c_str())) - return solids.Get(name.c_str()); + if (solids.Used(name)) + return solids[name]; else return NULL; } @@ -637,15 +637,15 @@ namespace netgen const SplineGeometry<2> * CSGeometry :: GetSplineCurve2d (const string & name) const { - if (splinecurves2d.Used(name.c_str())) - return splinecurves2d.Get(name.c_str()); + if (splinecurves2d.Used(name)) + return splinecurves2d[name]; else return NULL; } const SplineGeometry<3> * CSGeometry :: GetSplineCurve3d (const string & name) const { - if (splinecurves3d.Used(name.c_str())) - return splinecurves3d.Get(name.c_str()); + if (splinecurves3d.Used(name)) + return splinecurves3d[name]; else return NULL; } @@ -721,7 +721,7 @@ namespace netgen void CSGeometry :: SetFlags (const char * solidname, const Flags & flags) { - Solid * solid = solids.Elem(solidname); + Solid * solid = solids[solidname]; Array surfind; int i; diff --git a/libsrc/csg/csgeom.hpp b/libsrc/csg/csgeom.hpp index f82407b4..40d90d51 100644 --- a/libsrc/csg/csgeom.hpp +++ b/libsrc/csg/csgeom.hpp @@ -102,7 +102,7 @@ namespace netgen { private: /// all surfaces - SYMBOLTABLE surfaces; + SymbolTable surfaces; public: /// primitive of surface @@ -112,12 +112,12 @@ namespace netgen Array delete_them; /// all named solids - SYMBOLTABLE solids; + SymbolTable solids; /// all 2d splinecurves - SYMBOLTABLE< SplineGeometry<2>* > splinecurves2d; + SymbolTable< SplineGeometry<2>* > splinecurves2d; /// all 3d splinecurves - SYMBOLTABLE< SplineGeometry<3>* > splinecurves3d; + SymbolTable< SplineGeometry<3>* > splinecurves3d; /// all top level objects: solids and surfaces Array toplevelobjects; @@ -203,7 +203,7 @@ namespace netgen const Solid * GetSolid (const string & name) const; int GetNSolids () const { return solids.Size(); } const Solid * GetSolid (int i) const { return solids[i]; } - const SYMBOLTABLE & GetSolids () const { return solids; } + const SymbolTable & GetSolids () const { return solids; } void SetSplineCurve (const char * name, SplineGeometry<2> * spl); diff --git a/libsrc/csg/python_csg.cpp b/libsrc/csg/python_csg.cpp index 0db96e6d..c82db668 100644 --- a/libsrc/csg/python_csg.cpp +++ b/libsrc/csg/python_csg.cpp @@ -469,6 +469,17 @@ However, when r = 0, the top part becomes a point(tip) and meshing fails! self.AddSplineSurface(surf); }), py::arg("SplineSurface")) + .def("SingularFace", [] (CSGeometry & self, shared_ptr sol, shared_ptr surfaces, double factor) + { + int tlonum = -1; + for (int i = 0; i < self.GetNTopLevelObjects(); i++) + if (self.GetTopLevelObject(i)->GetSolid() == sol->GetSolid()) + tlonum = i; + if (tlonum == -1) throw NgException("not a top-level-object"); + if (!surfaces) surfaces = sol; + auto singface = new SingularFace(tlonum+1, surfaces->GetSolid(), factor); + self.singfaces.Append(singface); + }, py::arg("solid"), py::arg("surfaces")=nullptr, py::arg("factor")=0.25) .def("SingularEdge", [] (CSGeometry & self, shared_ptr s1,shared_ptr s2, double factor) { auto singedge = new SingularEdge(1, -1, self, s1->GetSolid(), s2->GetSolid(), factor); diff --git a/libsrc/csg/solid.cpp b/libsrc/csg/solid.cpp index a7e5fff5..c61f4e39 100644 --- a/libsrc/csg/solid.cpp +++ b/libsrc/csg/solid.cpp @@ -419,9 +419,9 @@ namespace netgen - static Solid * CreateSolidExpr (istream & ist, const SYMBOLTABLE & solids); - static Solid * CreateSolidTerm (istream & ist, const SYMBOLTABLE & solids); - static Solid * CreateSolidPrim (istream & ist, const SYMBOLTABLE & solids); + static Solid * CreateSolidExpr (istream & ist, const SymbolTable & solids); + static Solid * CreateSolidTerm (istream & ist, const SymbolTable & solids); + static Solid * CreateSolidPrim (istream & ist, const SymbolTable & solids); static void ReadString (istream & ist, char * str) { @@ -461,7 +461,7 @@ namespace netgen } - Solid * CreateSolidExpr (istream & ist, const SYMBOLTABLE & solids) + Solid * CreateSolidExpr (istream & ist, const SymbolTable & solids) { // cout << "create expr" << endl; @@ -484,7 +484,7 @@ namespace netgen return s1; } - Solid * CreateSolidTerm (istream & ist, const SYMBOLTABLE & solids) + Solid * CreateSolidTerm (istream & ist, const SymbolTable & solids) { // cout << "create term" << endl; @@ -508,7 +508,7 @@ namespace netgen return s1; } - Solid * CreateSolidPrim (istream & ist, const SYMBOLTABLE & solids) + Solid * CreateSolidPrim (istream & ist, const SymbolTable & solids) { Solid * s1; char ch; @@ -533,7 +533,7 @@ namespace netgen } (*testout) << "get terminal " << str << endl; - s1 = solids.Get(str); + s1 = solids[str]; if (s1) { // cout << "primitive: " << str << endl; @@ -545,7 +545,7 @@ namespace netgen } - Solid * Solid :: CreateSolid (istream & ist, const SYMBOLTABLE & solids) + Solid * Solid :: CreateSolid (istream & ist, const SymbolTable & solids) { Solid * nsol = CreateSolidExpr (ist, solids); nsol = new Solid (ROOT, nsol); diff --git a/libsrc/csg/solid.hpp b/libsrc/csg/solid.hpp index 4d620c04..193cb363 100644 --- a/libsrc/csg/solid.hpp +++ b/libsrc/csg/solid.hpp @@ -158,7 +158,7 @@ namespace netgen { return maxh; } void GetSolidData (ostream & ost, int first = 1) const; - static Solid * CreateSolid (istream & ist, const SYMBOLTABLE & solids); + static Solid * CreateSolid (istream & ist, const SymbolTable & solids); static BlockAllocator ball; diff --git a/libsrc/general/CMakeLists.txt b/libsrc/general/CMakeLists.txt index 00c59ea2..3f6c2aad 100644 --- a/libsrc/general/CMakeLists.txt +++ b/libsrc/general/CMakeLists.txt @@ -1,21 +1,19 @@ add_definitions(-DNGINTERFACE_EXPORTS) -add_library(gen OBJECT - array.cpp bitarray.cpp dynamicmem.cpp flags.cpp - hashtabl.cpp mystring.cpp ngexception.cpp optmem.cpp parthreads.cpp - profiler.cpp seti.cpp sort.cpp spbita2d.cpp symbolta.cpp table.cpp - mpi_interface.cpp gzstream.cpp +add_library(gen INTERFACE) +set(sdir ${CMAKE_CURRENT_SOURCE_DIR}) +target_sources(gen INTERFACE + ${sdir}/array.cpp ${sdir}/bitarray.cpp ${sdir}/dynamicmem.cpp ${sdir}/flags.cpp + ${sdir}/hashtabl.cpp ${sdir}/mystring.cpp ${sdir}/optmem.cpp ${sdir}/parthreads.cpp + ${sdir}/seti.cpp ${sdir}/sort.cpp ${sdir}/spbita2d.cpp ${sdir}/table.cpp + ${sdir}/mpi_interface.cpp ${sdir}/gzstream.cpp ) -set_target_properties( gen PROPERTIES POSITION_INDEPENDENT_CODE ON ) - -install( FILES ngexception.hpp DESTINATION ${NG_INSTALL_DIR_INCLUDE} COMPONENT netgen_devel ) - install(FILES array.hpp autodiff.hpp autoptr.hpp bitarray.hpp dynamicmem.hpp flags.hpp hashtabl.hpp mpi_interface.hpp myadt.hpp - ngsimd.hpp mystring.hpp netgenout.hpp ngexception.hpp ngpython.hpp - optmem.hpp parthreads.hpp profiler.hpp seti.hpp sort.hpp - spbita2d.hpp stack.hpp symbolta.hpp table.hpp template.hpp + ngsimd.hpp mystring.hpp netgenout.hpp ngpython.hpp + optmem.hpp parthreads.hpp seti.hpp sort.hpp + spbita2d.hpp stack.hpp table.hpp template.hpp gzstream.h DESTINATION ${NG_INSTALL_DIR_INCLUDE}/general COMPONENT netgen_devel ) diff --git a/libsrc/general/array.hpp b/libsrc/general/array.hpp index 16bb6e78..df0cd0a5 100644 --- a/libsrc/general/array.hpp +++ b/libsrc/general/array.hpp @@ -243,10 +243,11 @@ namespace netgen } explicit Array(size_t asize) - : FlatArray (asize, new T[asize]) + : FlatArray (asize, asize ? new T[asize] : nullptr) { - allocsize = asize; - ownmem = 1; + allocsize = asize; + if(asize) + ownmem = 1; } /// Generate array in user data diff --git a/libsrc/general/flags.cpp b/libsrc/general/flags.cpp index e5916f87..993413ed 100644 --- a/libsrc/general/flags.cpp +++ b/libsrc/general/flags.cpp @@ -83,7 +83,7 @@ namespace netgen Flags :: GetStringFlag (const char * name, const char * def) const { if (strflags.Used (name)) - return strflags.Get(name); + return strflags[name]; else return def; } @@ -91,7 +91,7 @@ namespace netgen double Flags :: GetNumFlag (const char * name, double def) const { if (numflags.Used (name)) - return numflags.Get(name); + return numflags[name]; else return def; } @@ -99,7 +99,7 @@ namespace netgen const double * Flags :: GetNumFlagPtr (const char * name) const { if (numflags.Used (name)) - return & ((SYMBOLTABLE&)numflags).Elem(name); + return & ((SymbolTable&)numflags)[name]; else return NULL; } @@ -107,7 +107,7 @@ namespace netgen double * Flags :: GetNumFlagPtr (const char * name) { if (numflags.Used (name)) - return & ((SYMBOLTABLE&)numflags).Elem(name); + return & ((SymbolTable&)numflags)[name]; else return NULL; } @@ -122,7 +122,7 @@ namespace netgen Flags :: GetStringListFlag (const char * name) const { if (strlistflags.Used (name)) - return *strlistflags.Get(name); + return *strlistflags[name]; else { static Array dummy_array(0); @@ -134,7 +134,7 @@ namespace netgen Flags ::GetNumListFlag (const char * name) const { if (numlistflags.Used (name)) - return *numlistflags.Get(name); + return *numlistflags[name]; else { static Array dummy_array(0); @@ -170,9 +170,9 @@ namespace netgen ofstream outfile (filename); for (i = 1; i <= strflags.Size(); i++) - outfile << strflags.GetName(i) << " = " << strflags.Get(i) << endl; + outfile << strflags.GetName(i) << " = " << strflags[i] << endl; for (i = 1; i <= numflags.Size(); i++) - outfile << numflags.GetName(i) << " = " << numflags.Get(i) << endl; + outfile << numflags.GetName(i) << " = " << numflags[i] << endl; for (i = 1; i <= defflags.Size(); i++) outfile << defflags.GetName(i) << endl; } @@ -184,9 +184,9 @@ namespace netgen int i; for (i = 1; i <= strflags.Size(); i++) - ost << strflags.GetName(i) << " = " << strflags.Get(i) << endl; + ost << strflags.GetName(i) << " = " << strflags[i] << endl; for (i = 1; i <= numflags.Size(); i++) - ost << numflags.GetName(i) << " = " << numflags.Get(i) << endl; + ost << numflags.GetName(i) << " = " << numflags[i] << endl; for (i = 1; i <= defflags.Size(); i++) ost << defflags.GetName(i) << endl; } diff --git a/libsrc/general/flags.hpp b/libsrc/general/flags.hpp index e156d328..e7ba0382 100644 --- a/libsrc/general/flags.hpp +++ b/libsrc/general/flags.hpp @@ -19,15 +19,15 @@ namespace netgen class Flags { /// - SYMBOLTABLE strflags; + SymbolTable strflags; /// - SYMBOLTABLE numflags; + SymbolTable numflags; /// - SYMBOLTABLE defflags; + SymbolTable defflags; /// - SYMBOLTABLE*> strlistflags; + SymbolTable*> strlistflags; /// - SYMBOLTABLE*> numlistflags; + SymbolTable*> numlistflags; public: /// DLL_HEADER Flags (); diff --git a/libsrc/general/mpi_interface.hpp b/libsrc/general/mpi_interface.hpp index c45494b6..eb79865e 100644 --- a/libsrc/general/mpi_interface.hpp +++ b/libsrc/general/mpi_interface.hpp @@ -14,101 +14,172 @@ namespace netgen { + using ngcore::id; + using ngcore::ntasks; - extern DLL_HEADER int id, ntasks; - +#ifndef PARALLEL + /** without MPI, we need a dummy typedef **/ + typedef int MPI_Comm; +#endif + + /** This is the "standard" communicator that will be used for netgen-objects. **/ + extern DLL_HEADER MPI_Comm ng_comm; + +#ifdef PARALLEL + inline int MyMPI_GetNTasks (MPI_Comm comm = ng_comm) + { + int ntasks; + MPI_Comm_size(comm, &ntasks); + return ntasks; + } + inline int MyMPI_GetId (MPI_Comm comm = ng_comm) + { + int id; + MPI_Comm_rank(comm, &id); + return id; + } +#else + enum { MPI_COMM_WORLD = 12345, MPI_COMM_NULL = 0}; + inline int MyMPI_GetNTasks (MPI_Comm comm = ng_comm) { return 1; } + inline int MyMPI_GetId (MPI_Comm comm = ng_comm) { return 0; } +#endif + +#ifdef PARALLEL + // For python wrapping of communicators + struct PyMPI_Comm { + MPI_Comm comm; + bool owns_comm; + PyMPI_Comm (MPI_Comm _comm, bool _owns_comm = false) : comm(_comm), owns_comm(_owns_comm) { } + PyMPI_Comm (const PyMPI_Comm & c) = delete; + ~PyMPI_Comm () { + if (owns_comm) + MPI_Comm_free(&comm); + } + inline int Rank() const { return MyMPI_GetId(comm); } + inline int Size() const { return MyMPI_GetNTasks(comm); } + }; +#else + // dummy without MPI + struct PyMPI_Comm { + MPI_Comm comm = 0; + PyMPI_Comm (MPI_Comm _comm, bool _owns_comm = false) { } + ~PyMPI_Comm () { } + inline int Rank() const { return 0; } + inline int Size() const { return 1; } + }; +#endif + +#ifdef PARALLEL + template + inline MPI_Datatype MyGetMPIType ( ) + { cerr << "ERROR in GetMPIType() -- no type found" << endl;return 0; } + template <> + inline MPI_Datatype MyGetMPIType ( ) + { return MPI_INT; } + template <> + inline MPI_Datatype MyGetMPIType ( ) + { return MPI_DOUBLE; } + template <> + inline MPI_Datatype MyGetMPIType ( ) + { return MPI_CHAR; } + template<> + inline MPI_Datatype MyGetMPIType ( ) + { return MPI_UINT64_T; } +#else + typedef int MPI_Datatype; + template inline MPI_Datatype MyGetMPIType ( ) { return 0; } +#endif + +#ifdef PARALLEL + inline MPI_Comm MyMPI_SubCommunicator(MPI_Comm comm, Array & procs) + { + MPI_Comm subcomm; + MPI_Group gcomm, gsubcomm; + MPI_Comm_group(comm, &gcomm); + MPI_Group_incl(gcomm, procs.Size(), &(procs[0]), &gsubcomm); + MPI_Comm_create_group(comm, gsubcomm, 6969, &subcomm); + return subcomm; + } +#else + inline MPI_Comm MyMPI_SubCommunicator(MPI_Comm comm, Array & procs) + { return comm; } +#endif #ifdef PARALLEL - enum { MPI_TAG_CMD = 110 }; enum { MPI_TAG_MESH = 210 }; enum { MPI_TAG_VIS = 310 }; - extern MPI_Comm mesh_comm; - - template - MPI_Datatype MyGetMPIType ( ) - { cerr << "ERROR in GetMPIType() -- no type found" << endl;return 0; } - - template <> - inline MPI_Datatype MyGetMPIType ( ) - { return MPI_INT; } - - template <> - inline MPI_Datatype MyGetMPIType ( ) - { return MPI_DOUBLE; } - - - inline void MyMPI_Send (int i, int dest, int tag) + inline void MyMPI_Send (int i, int dest, int tag, MPI_Comm comm = ng_comm) { int hi = i; - MPI_Send( &hi, 1, MPI_INT, dest, tag, MPI_COMM_WORLD); + MPI_Send( &hi, 1, MPI_INT, dest, tag, comm); } - inline void MyMPI_Recv (int & i, int src, int tag) + inline void MyMPI_Recv (int & i, int src, int tag, MPI_Comm comm = ng_comm) { MPI_Status status; - MPI_Recv( &i, 1, MPI_INT, src, tag, MPI_COMM_WORLD, &status); + MPI_Recv( &i, 1, MPI_INT, src, tag, comm, &status); } - inline void MyMPI_Send (const string & s, int dest, int tag) + inline void MyMPI_Send (const string & s, int dest, int tag, MPI_Comm comm = ng_comm) { - MPI_Send( const_cast (s.c_str()), s.length(), MPI_CHAR, dest, tag, MPI_COMM_WORLD); + MPI_Send( const_cast (s.c_str()), s.length(), MPI_CHAR, dest, tag, comm); } - inline void MyMPI_Recv (string & s, int src, int tag) + inline void MyMPI_Recv (string & s, int src, int tag, MPI_Comm comm = ng_comm) { MPI_Status status; int len; MPI_Probe (src, tag, MPI_COMM_WORLD, &status); MPI_Get_count (&status, MPI_CHAR, &len); s.assign (len, ' '); - MPI_Recv( &s[0], len, MPI_CHAR, src, tag, MPI_COMM_WORLD, &status); + MPI_Recv( &s[0], len, MPI_CHAR, src, tag, comm, &status); } template - inline void MyMPI_Send (FlatArray s, int dest, int tag) + inline void MyMPI_Send (FlatArray s, int dest, int tag, MPI_Comm comm = ng_comm) { - MPI_Send( &s.First(), s.Size(), MyGetMPIType(), dest, tag, MPI_COMM_WORLD); + MPI_Send( &s.First(), s.Size(), MyGetMPIType(), dest, tag, comm); } template - inline void MyMPI_Recv ( FlatArray s, int src, int tag) + inline void MyMPI_Recv ( FlatArray s, int src, int tag, MPI_Comm comm = ng_comm) { MPI_Status status; - MPI_Recv( &s.First(), s.Size(), MyGetMPIType(), src, tag, MPI_COMM_WORLD, &status); + MPI_Recv( &s.First(), s.Size(), MyGetMPIType(), src, tag, comm, &status); } template - inline void MyMPI_Recv ( Array & s, int src, int tag) + inline void MyMPI_Recv ( Array & s, int src, int tag, MPI_Comm comm = ng_comm) { MPI_Status status; int len; - MPI_Probe (src, tag, MPI_COMM_WORLD, &status); + MPI_Probe (src, tag, comm, &status); MPI_Get_count (&status, MyGetMPIType(), &len); s.SetSize (len); - MPI_Recv( &s.First(), len, MyGetMPIType(), src, tag, MPI_COMM_WORLD, &status); + MPI_Recv( &s.First(), len, MyGetMPIType(), src, tag, comm, &status); } template - inline int MyMPI_Recv ( Array & s, int tag) + inline int MyMPI_Recv ( Array & s, int tag, MPI_Comm comm = ng_comm) { MPI_Status status; int len; - MPI_Probe (MPI_ANY_SOURCE, tag, MPI_COMM_WORLD, &status); + MPI_Probe (MPI_ANY_SOURCE, tag, comm, &status); int src = status.MPI_SOURCE; MPI_Get_count (&status, MyGetMPIType(), &len); s.SetSize (len); - MPI_Recv( &s.First(), len, MyGetMPIType(), src, tag, MPI_COMM_WORLD, &status); + MPI_Recv( &s.First(), len, MyGetMPIType(), src, tag, comm, &status); return src; } @@ -130,7 +201,7 @@ namespace netgen */ template - inline MPI_Request MyMPI_ISend (FlatArray s, int dest, int tag, MPI_Comm comm = MPI_COMM_WORLD) + inline MPI_Request MyMPI_ISend (FlatArray s, int dest, int tag, MPI_Comm comm = ng_comm) { MPI_Request request; MPI_Isend( &s.First(), s.Size(), MyGetMPIType(), dest, tag, comm, &request); @@ -139,7 +210,7 @@ namespace netgen template - inline MPI_Request MyMPI_IRecv (FlatArray s, int dest, int tag, MPI_Comm comm = MPI_COMM_WORLD) + inline MPI_Request MyMPI_IRecv (FlatArray s, int dest, int tag, MPI_Comm comm = ng_comm) { MPI_Request request; MPI_Irecv( &s.First(), s.Size(), MyGetMPIType(), dest, tag, comm, &request); @@ -204,11 +275,10 @@ namespace netgen template inline void MyMPI_ExchangeTable (TABLE & send_data, TABLE & recv_data, int tag, - MPI_Comm comm = MPI_COMM_WORLD) + MPI_Comm comm = ng_comm) { - int ntasks, rank; - MPI_Comm_size(comm, &ntasks); - MPI_Comm_rank(comm, &rank); + int rank = MyMPI_GetId(comm); + int ntasks = MyMPI_GetNTasks(comm); Array send_sizes(ntasks); Array recv_sizes(ntasks); @@ -252,22 +322,22 @@ namespace netgen template - inline void MyMPI_Bcast (T & s, MPI_Comm comm = MPI_COMM_WORLD) + inline void MyMPI_Bcast (T & s, MPI_Comm comm = ng_comm) { MPI_Bcast (&s, 1, MyGetMPIType(), 0, comm); } template - inline void MyMPI_Bcast (Array & s, MPI_Comm comm = MPI_COMM_WORLD) + inline void MyMPI_Bcast (Array & s, MPI_Comm comm = ng_comm) { int size = s.Size(); MyMPI_Bcast (size, comm); - if (id != 0) s.SetSize (size); + if (MyMPI_GetId(comm) != 0) s.SetSize (size); MPI_Bcast (&s[0], size, MyGetMPIType(), 0, comm); } template - inline void MyMPI_Bcast (Array & s, int root, MPI_Comm comm = MPI_COMM_WORLD) + inline void MyMPI_Bcast (Array & s, int root, MPI_Comm comm = ng_comm) { int id; MPI_Comm_rank(comm, &id); @@ -280,19 +350,19 @@ namespace netgen } template - inline void MyMPI_Allgather (const T & send, FlatArray recv, MPI_Comm comm) + inline void MyMPI_Allgather (const T & send, FlatArray recv, MPI_Comm comm = ng_comm) { MPI_Allgather( const_cast (&send), 1, MyGetMPIType(), &recv[0], 1, MyGetMPIType(), comm); } template - inline void MyMPI_Alltoall (FlatArray send, FlatArray recv, MPI_Comm comm) + inline void MyMPI_Alltoall (FlatArray send, FlatArray recv, MPI_Comm comm = ng_comm) { MPI_Alltoall( &send[0], 1, MyGetMPIType(), &recv[0], 1, MyGetMPIType(), comm); } // template -// inline void MyMPI_Alltoall_Block (FlatArray send, FlatArray recv, int blocklen, MPI_Comm comm) +// inline void MyMPI_Alltoall_Block (FlatArray send, FlatArray recv, int blocklen, MPI_Comm comm = ng_comm) // { // MPI_Alltoall( &send[0], blocklen, MyGetMPIType(), &recv[0], blocklen, MyGetMPIType(), comm); // } diff --git a/libsrc/general/myadt.hpp b/libsrc/general/myadt.hpp index 601a3da0..4284a42f 100644 --- a/libsrc/general/myadt.hpp +++ b/libsrc/general/myadt.hpp @@ -21,8 +21,8 @@ namespace netgen { using namespace ngcore; + using NgException = Exception; } -#include "ngexception.hpp" #include "parthreads.hpp" // #include "moveablemem.hpp" #include "dynamicmem.hpp" @@ -33,7 +33,6 @@ namespace netgen #include "hashtabl.hpp" -#include "symbolta.hpp" #include "bitarray.hpp" #include "flags.hpp" #include "spbita2d.hpp" @@ -44,7 +43,6 @@ namespace netgen #include "sort.hpp" #include "stack.hpp" #include "mystring.hpp" -#include "profiler.hpp" #include "mpi_interface.hpp" #include "netgenout.hpp" diff --git a/libsrc/general/netgenout.hpp b/libsrc/general/netgenout.hpp index fafe495e..076a6df1 100644 --- a/libsrc/general/netgenout.hpp +++ b/libsrc/general/netgenout.hpp @@ -4,14 +4,11 @@ // #include // #include // #include +#include "mpi_interface.hpp" namespace netgen { -#ifdef PARALLEL -extern int id; -extern int ntasks; -#endif DLL_HEADER extern int printmessage_importance; DLL_HEADER extern int printdots; diff --git a/libsrc/general/ngexception.cpp b/libsrc/general/ngexception.cpp deleted file mode 100644 index 6bcd2cc9..00000000 --- a/libsrc/general/ngexception.cpp +++ /dev/null @@ -1,33 +0,0 @@ -/**************************************************************************/ -/* File: ngexception.cpp */ -/* Author: Joachim Schoeberl */ -/* Date: 16. Jan. 02 */ -/**************************************************************************/ - -#include - -namespace netgen -{ - //using namespace netgen; - - - - NgException :: NgException (const string & s) - : m_what(s) - { - ; - } - - - NgException :: ~NgException () - { - ; - } - - /// append string to description - void NgException :: Append (const string & s) - { - m_what += s; - } - -} diff --git a/libsrc/general/ngexception.hpp b/libsrc/general/ngexception.hpp deleted file mode 100644 index 6e06498c..00000000 --- a/libsrc/general/ngexception.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef FILE_NGEXCEPTION -#define FILE_NGEXCEPTION - -/**************************************************************************/ -/* File: ngexception.hpp */ -/* Author: Joachim Schoeberl */ -/* Date: 16. Jan. 2002 */ -/**************************************************************************/ - -namespace netgen -{ - -/// Base class for all ng exceptions -class NgException : public std::exception -{ - /// verbal description of exception - string m_what; -public: - /// - DLL_HEADER NgException (const string & s); - /// - DLL_HEADER virtual ~NgException (); - - /// append string to description - DLL_HEADER void Append (const string & s); - // void Append (const char * s); - - /// verbal description of exception - const string & What() const { return m_what; } - virtual const char* what() const noexcept override { return m_what.c_str(); } -}; -} - -#endif diff --git a/libsrc/general/ngpython.hpp b/libsrc/general/ngpython.hpp index fa9862b1..722709c5 100644 --- a/libsrc/general/ngpython.hpp +++ b/libsrc/general/ngpython.hpp @@ -3,6 +3,7 @@ #include #include #include +#include namespace py = pybind11; #include #include @@ -66,16 +67,7 @@ namespace netgen return static_cast::pointer>(lambda); } - - template - inline std::string ToString (const T& t) - { - std::stringstream ss; - ss << t; - return ss.str(); - } - -} +} // namespace netgen #endif diff --git a/libsrc/general/profiler.cpp b/libsrc/general/profiler.cpp deleted file mode 100644 index e6a9685b..00000000 --- a/libsrc/general/profiler.cpp +++ /dev/null @@ -1,131 +0,0 @@ -/**************************************************************************/ -/* File: profiler.cpp */ -/* Author: Joachim Schoeberl */ -/* Date: 19. Apr. 2002 */ -/**************************************************************************/ - - -#include - -namespace netgen -{ - //using namespace netgen; - - long int NgProfiler::tottimes[SIZE]; - long int NgProfiler::starttimes[SIZE]; - long int NgProfiler::counts[SIZE]; - string NgProfiler::names[SIZE]; - int NgProfiler::usedcounter[SIZE]; - - - NgProfiler :: NgProfiler() - { - for (int i = 0; i < SIZE; i++) - { - tottimes[i] = 0; - usedcounter[i] = 0; - } - - total_timer = CreateTimer ("total CPU time"); - StartTimer (total_timer); - } - - NgProfiler :: ~NgProfiler() - { -#ifndef PARALLEL - StopTimer (total_timer); -#endif - - //ofstream prof; - //prof.open("ng.prof"); - - // ofstream-constructor may be called after STL-stuff is destructed, - // which leads to an "order of destruction"-problem, - // thus we use the C-variant: - - if (getenv ("NGPROFILE")) - { - char filename[100]; -#ifdef PARALLEL - sprintf (filename, "netgen.prof.%d", id); -#else - sprintf (filename, "netgen.prof"); -#endif - - if (id == 0) printf ("write profile to file netgen.prof\n"); - FILE *prof = fopen(filename,"w"); - Print (prof); - fclose(prof); - } - } - - -// void NgProfiler :: Print (ostream & prof) -// { -// for (int i = 0; i < SIZE; i++) -// if (counts[i] != 0 || usedcounter[i] != 0) -// { -// prof.setf (ios::fixed, ios::floatfield); -// prof.setf (ios::showpoint); - -// prof // << "job " << setw(3) << i -// << "calls " << setw(8) << counts[i] -// << ", time " << setprecision(2) << setw(6) << double(tottimes[i]) / CLOCKS_PER_SEC << " sec"; - -// if (usedcounter[i]) -// prof << " " << names[i]; -// else -// prof << " " << i; - -// prof << endl; -// } -// } - - - void NgProfiler :: Print (FILE * prof) - { - for (int i = 0; i < SIZE; i++) - if (counts[i] != 0 || usedcounter[i] != 0) - { - //fprintf(prof,"job %3i calls %8i, time %6.2f sec",i,counts[i],double(tottimes[i]) / CLOCKS_PER_SEC); -#ifndef USE_TSC - fprintf(prof,"calls %8li, time %6.2f sec",counts[i],double(tottimes[i]) / CLOCKS_PER_SEC); -#else - fprintf(prof,"calls %8li, time %6.2f sec",counts[i],double(tottimes[i]) / 2.7e9); -#endif - if(usedcounter[i]) - fprintf(prof," %s",names[i].c_str()); - else - fprintf(prof," %i",i); - fprintf(prof,"\n"); - } - } - - int NgProfiler :: CreateTimer (const string & name) - { - for (int i = SIZE-1; i > 0; i--) - if(names[i] == name) - return i; - - for (int i = SIZE-1; i > 0; i--) - if (!usedcounter[i]) - { - usedcounter[i] = 1; - names[i] = name; - return i; - } - return -1; - } - - - void NgProfiler :: ClearTimers () - { - for (int i = 0; i < SIZE; i++) - { - tottimes[i] = 0; - counts[i] = 0; - } - } - - NgProfiler prof; -} diff --git a/libsrc/general/profiler.hpp b/libsrc/general/profiler.hpp deleted file mode 100644 index 98039114..00000000 --- a/libsrc/general/profiler.hpp +++ /dev/null @@ -1,88 +0,0 @@ -#ifndef FILE_NG_PROFILER -#define FILE_NG_PROFILER - -/**************************************************************************/ -/* File: profiler.hpp */ -/* Author: Joachim Schoeberl */ -/* Date: 5. Jan. 2005 */ -/**************************************************************************/ - - - -#ifdef VTRACE -#include "vt_user.h" -#else - #define VT_USER_START(n) - #define VT_USER_END(n) - #define VT_TRACER(n) -#endif - - -// #define USE_TSC -#ifdef USE_TSC -#include // for __rdtsc() CPU time step counter -#endif - -namespace netgen -{ - -class NgProfiler -{ - enum { SIZE = 1000 }; - - static long int tottimes[SIZE]; - static long int starttimes[SIZE]; - static long int counts[SIZE]; - static string names[SIZE]; - static int usedcounter[SIZE]; - - int total_timer; -public: - NgProfiler(); - ~NgProfiler(); - static int CreateTimer (const string & name); -#ifndef USE_TSC - static void StartTimer (int nr) - { - starttimes[nr] = clock(); counts[nr]++; - // VT_USER_START (const_cast (names[nr].c_str())); - VT_USER_START ( (char * const) (names[nr].c_str())); - } - static void StopTimer (int nr) - { - tottimes[nr] += clock()-starttimes[nr]; - VT_USER_END (const_cast (names[nr].c_str())); - } -#else - static void StartTimer (int nr) - { - starttimes[nr] = __rdtsc(); counts[nr]++; - // VT_USER_START (const_cast (names[nr].c_str())); - // VT_USER_START ( (char * const) (names[nr].c_str())); - } - static void StopTimer (int nr) - { - tottimes[nr] += __rdtsc()-starttimes[nr]; - VT_USER_END (const_cast (names[nr].c_str())); - } -#endif - - - //static void Print (ostream & ost); - static void Print (FILE * prof); - - static void ClearTimers (); - - class RegionTimer - { - int nr; - public: - RegionTimer (int anr) : nr(anr) - { StartTimer (nr); } - ~RegionTimer () { StopTimer (nr); } - }; -}; - -} - -#endif diff --git a/libsrc/general/symbolta.cpp b/libsrc/general/symbolta.cpp deleted file mode 100644 index bd35ac7c..00000000 --- a/libsrc/general/symbolta.cpp +++ /dev/null @@ -1,52 +0,0 @@ -/**************************************************************************/ -/* File: symbolta.cc */ -/* Author: Joachim Schoeberl */ -/* Date: 01. Jun. 95 */ -/**************************************************************************/ - -/* - Abstract data type Symbol Table -*/ - -#include -#include - - -#ifndef FILE_SYMBOLTABLECC -#define FILE_SYMBOLTABLECC -// necessary for SGI ???? - - -namespace netgen -{ - //using namespace netgen; - - BASE_SYMBOLTABLE :: BASE_SYMBOLTABLE () - { - ; - } - - - BASE_SYMBOLTABLE :: ~BASE_SYMBOLTABLE() - { - DelNames(); - } - - - void BASE_SYMBOLTABLE :: DelNames() - { - for (int i = 0; i < names.Size(); i++) - delete [] names[i]; - names.SetSize (0); - } - - int BASE_SYMBOLTABLE :: Index (const char * name) const - { - if (!name) return 0; - for (int i = 0; i < names.Size(); i++) - if (strcmp (names[i], name) == 0) return i+1; - return 0; - } -} - -#endif diff --git a/libsrc/general/symbolta.hpp b/libsrc/general/symbolta.hpp deleted file mode 100644 index b599ea42..00000000 --- a/libsrc/general/symbolta.hpp +++ /dev/null @@ -1,163 +0,0 @@ -#ifndef FILE_SYMBOLTA -#define FILE_SYMBOLTA - - -/**************************************************************************/ -/* File: symbolta.hh */ -/* Author: Joachim Schoeberl */ -/* Date: 01. Jun. 95 */ -/**************************************************************************/ - -namespace netgen -{ - -/** - Base class for the generic SYMBOLTABLE. - An array of identifiers is maintained. -*/ -class DLL_HEADER BASE_SYMBOLTABLE -{ -protected: - /// identifiers - Array names; - -public: - /// Constructor - BASE_SYMBOLTABLE (); - /// - ~BASE_SYMBOLTABLE (); - /// - void DelNames (); - /// Index of symbol name, returns 0 if not used. - int Index (const char * name) const; -}; - - -/** - Abstract data type Symbol Table. - - To a string an value of the generic type T is associated. - The string is not copied into the symbol table class! -*/ -template -class SYMBOLTABLE : public BASE_SYMBOLTABLE -{ -private: - /// Associated data - Array data; - -public: - /// Creates a symboltable - inline SYMBOLTABLE (); - /// Returns size of symboltable - inline INDEX Size() const; - /// Returns reference to element, error if not used - inline T & Elem (const char * name); - /// Returns reference to i-th element - inline T & Elem (int i) - { return data.Elem(i); } - /// Returns element, error if not used - inline const T & Get (const char * name) const; - /// Returns i-th element - inline const T & Get (int i) const; - /// Returns name of i-th element - inline const char* GetName (int i) const; - /// Associates el to the string name, overrides if name is used - inline void Set (const char * name, const T & el); - /// Checks whether name is used - inline bool Used (const char * name) const; - /// Deletes symboltable - inline void DeleteAll (); - - void DoArchive(Archive& archive) { archive & names & data;} - - inline T & operator[] (int i) - { return data[i]; } - inline const T & operator[] (int i) const - { return data[i]; } - -private: - /// Prevents from copying symboltable by pointer assignment - SYMBOLTABLE & operator= (SYMBOLTABLE &); -}; - - - - -template -inline SYMBOLTABLE :: SYMBOLTABLE () -{ - ; -} - - -template -inline INDEX SYMBOLTABLE :: Size() const -{ - return data.Size(); -} - -template -inline T & SYMBOLTABLE :: Elem (const char * name) -{ - int i = Index (name); - if (i) - return data.Elem (i); - else - return data.Elem(1); -} - -template -inline const T & SYMBOLTABLE :: Get (const char * name) const -{ - int i; - i = Index (name); - if (i) - return data.Get(i); - else - return data.Get(1); -} - -template -inline const T & SYMBOLTABLE :: Get (int i) const -{ - return data.Get(i); -} - -template -inline const char* SYMBOLTABLE :: GetName (int i) const -{ - return names.Get(i); -} - -template -inline void SYMBOLTABLE :: Set (const char * name, const T & el) -{ - int i; - i = Index (name); - if (i) - data.Set(i, el); - else - { - data.Append (el); - char * hname = new char [strlen (name) + 1]; - strcpy (hname, name); - names.Append (hname); - } -} - -template -inline bool SYMBOLTABLE :: Used (const char * name) const -{ - return (Index(name)) ? true : false; -} - -template -inline void SYMBOLTABLE :: DeleteAll () -{ - DelNames(); - data.DeleteAll(); -} - -} -#endif diff --git a/libsrc/geom2d/CMakeLists.txt b/libsrc/geom2d/CMakeLists.txt index 754c79fb..43c619a0 100644 --- a/libsrc/geom2d/CMakeLists.txt +++ b/libsrc/geom2d/CMakeLists.txt @@ -4,10 +4,10 @@ if(APPLE) set_target_properties( geom2d PROPERTIES SUFFIX ".so") endif(APPLE) -if(NOT WIN32) - target_link_libraries(geom2d mesh ${PYTHON_LIBRARIES}) - install( TARGETS geom2d ${NG_INSTALL_DIR}) -endif(NOT WIN32) +target_link_libraries(geom2d mesh ${PYTHON_LIBRARIES}) +install( TARGETS geom2d ${NG_INSTALL_DIR}) + +target_link_libraries(geom2d ngcore) if(USE_GUI) add_library(geom2dvis ${NG_LIB_TYPE} vsgeom2d.cpp) diff --git a/libsrc/geom2d/genmesh2d.cpp b/libsrc/geom2d/genmesh2d.cpp index 43193e3a..fde8b16c 100644 --- a/libsrc/geom2d/genmesh2d.cpp +++ b/libsrc/geom2d/genmesh2d.cpp @@ -55,7 +55,8 @@ namespace netgen sum += dt / fun; } - int nel = int (sum+1); + int nel = int (sum+0.5); + if (nel == 0) nel = 1; fperel = sum / nel; points.Append (0); @@ -457,9 +458,9 @@ namespace netgen for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++) if ( (*mesh)[si].si > maxsegmentindex) maxsegmentindex = (*mesh)[si].si; - mesh->SetNBCNames(maxsegmentindex+1); + mesh->SetNBCNames(maxsegmentindex); - for ( int sindex = 0; sindex <= maxsegmentindex; sindex++ ) + for ( int sindex = 0; sindex < maxsegmentindex; sindex++ ) mesh->SetBCName ( sindex, geometry.GetBCName( sindex+1 ) ); for (SegmentIndex si = 0; si < mesh->GetNSeg(); si++) diff --git a/libsrc/geom2d/geometry2d.hpp b/libsrc/geom2d/geometry2d.hpp index 52e7a20e..0d4bcb00 100644 --- a/libsrc/geom2d/geometry2d.hpp +++ b/libsrc/geom2d/geometry2d.hpp @@ -175,6 +175,7 @@ namespace netgen void CopyEdgeMesh (int from, int to, Mesh & mesh2d, Point3dTree & searchtree); + size_t GetNDomains() const { return materials.Size(); } void GetMaterial (int domnr, char* & material ); void SetMaterial (int domnr, const string & material); diff --git a/libsrc/geom2d/python_geom2d.cpp b/libsrc/geom2d/python_geom2d.cpp index a16c13f7..f2b20127 100644 --- a/libsrc/geom2d/python_geom2d.cpp +++ b/libsrc/geom2d/python_geom2d.cpp @@ -15,6 +15,22 @@ namespace netgen DLL_HEADER void ExportGeom2d(py::module &m) { + py::class_> + (m, "Spline", "Spline of a SplineGeometry object") + .def_property("leftdom", [] (SplineSegExt& self) { return self.leftdom; }, + [](SplineSegExt& self, int dom) { self.leftdom = dom; }) + .def_property("rightdom", [] (SplineSegExt& self) { return self.rightdom; }, + [](SplineSegExt& self, int dom) { self.rightdom = dom; }) + .def_property_readonly("bc", [] (SplineSegExt& self) { return self.bc; }) + .def("GetNormal", [](SplineSegExt& self, double t) + { + auto tang = self.GetTangent(t).Normalize(); + return Vec<2>(tang[1], -tang[0]); + }) + .def("StartPoint", [](SplineSegExt& self) { return Point<2>(self.StartPI()); }) + .def("EndPoint", [](SplineSegExt& self) { return Point<2>(self.EndPI()); }) + ; + py::class_> (m, "SplineGeometry", "a 2d boundary representation geometry model by lines and splines", @@ -121,38 +137,21 @@ DLL_HEADER void ExportGeom2d(py::module &m) seg->copyfrom = -1; self.AppendSegment(seg); }), py::arg("point_indices"), py::arg("leftdomain") = 1, py::arg("rightdomain") = py::int_(0)) - //.def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, int point_index1, int point_index2)//, int leftdomain, int rightdomain) - // { - // LineSeg<2> * l = new LineSeg<2>(self.GetPoint(point_index1), self.GetPoint(point_index2)); - // SplineSegExt * seg = new SplineSegExt(*l); - // seg->leftdom = 1;// leftdomain; - // seg->rightdom = 0;// rightdomain; - // seg->hmax = 1e99; - // seg->reffak = 1; - // seg->copyfrom = -1; - - // self.AppendSegment(seg); - // }))//, (py::arg("self"), py::arg("point_index1"), py::arg("point_index2"), py::arg("leftdomain") = 1, py::arg("rightdomain") = 0) ) - //.def("AppendSegment", FunctionPointer([](SplineGeometry2d &self, int point_index1, int point_index2, int point_index3)//, int leftdomain, int rightdomain) - // { - // SplineSeg3<2> * seg3 = new SplineSeg3<2>(self.GetPoint(point_index1), self.GetPoint(point_index2), self.GetPoint(point_index3)); - // SplineSegExt * seg = new SplineSegExt(*seg3); - // seg->leftdom = 1;// leftdomain; - // seg->rightdom = 0;// rightdomain; - // seg->hmax = 1e99; - // seg->reffak = 1; - // seg->copyfrom = -1; - // self.AppendSegment(seg); - // }))//, (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("SetDomainMaxH", &SplineGeometry2d::SetDomainMaxh) + .def("GetBCName", [](SplineGeometry2d& self, size_t index) { return self.GetBCName(index); }) + .def("GetNDomains", [](SplineGeometry2d& self) { return self.GetNDomains(); }) + .def("GetNSplines", [](SplineGeometry2d& self) { return self.splines.Size(); }) + .def("GetSpline", [](SplineGeometry2d& self, size_t index) + { return shared_ptr(&self.GetSpline(index), NOOP_Deleter); }, + py::return_value_policy::reference_internal) + .def("GetNPoints", [](SplineGeometry2d& self) { return self.GetNP(); }) + .def("GetPoint", [](SplineGeometry2d& self, size_t index) { return Point<2>(self.GetPoint(index)); }) - .def("PlotData", FunctionPointer([](SplineGeometry2d &self) { Box<2> box(self.GetBoundingBox()); diff --git a/libsrc/gprim/CMakeLists.txt b/libsrc/gprim/CMakeLists.txt index 3957476a..402973dc 100644 --- a/libsrc/gprim/CMakeLists.txt +++ b/libsrc/gprim/CMakeLists.txt @@ -1,11 +1,11 @@ add_definitions(-DNGINTERFACE_EXPORTS) -add_library(gprim OBJECT - adtree.cpp geom2d.cpp geom3d.cpp geomfuncs.cpp - geomtest3d.cpp transform3d.cpp spline.cpp splinegeometry.cpp +add_library(gprim INTERFACE) +set(sdir ${CMAKE_CURRENT_SOURCE_DIR}) +target_sources(gprim INTERFACE + ${sdir}/adtree.cpp ${sdir}/geom2d.cpp ${sdir}/geom3d.cpp ${sdir}/geomfuncs.cpp + ${sdir}/geomtest3d.cpp ${sdir}/transform3d.cpp ${sdir}/spline.cpp ${sdir}/splinegeometry.cpp ) -set_target_properties( gprim PROPERTIES POSITION_INDEPENDENT_CODE ON ) - install(FILES adtree.hpp geom2d.hpp geom3d.hpp geomfuncs.hpp geomobjects2.hpp geomobjects.hpp geomops2.hpp geomops.hpp geomtest3d.hpp gprim.hpp diff --git a/libsrc/gprim/adtree.cpp b/libsrc/gprim/adtree.cpp index 41237acc..cba845a7 100644 --- a/libsrc/gprim/adtree.cpp +++ b/libsrc/gprim/adtree.cpp @@ -2400,13 +2400,13 @@ namespace netgen Array & pis) const { Point<2*dim> tpmin, tpmax; - + double tol = Tolerance(); for (size_t i = 0; i < dim; i++) { tpmin(i) = boxpmin(i); - tpmax(i) = pmax(i); + tpmax(i) = pmax(i)+tol; - tpmin(i+dim) = pmin(i); + tpmin(i+dim) = pmin(i)-tol; tpmax(i+dim) = boxpmax(i); } diff --git a/libsrc/gprim/adtree.hpp b/libsrc/gprim/adtree.hpp index 174e9d9a..11694b67 100644 --- a/libsrc/gprim/adtree.hpp +++ b/libsrc/gprim/adtree.hpp @@ -570,9 +570,9 @@ public: { tree->DeleteElement(pi); } void GetIntersecting (const Point & pmin, const Point & pmax, Array & pis) const; - - // const T_ADTree<2*dim> & Tree() const { return *tree; }; - // T_ADTree<2*dim> & Tree() { return *tree; }; + double Tolerance() const { return 1e-7 * Dist(boxpmax, boxpmin); } // single precision + const auto & Tree() const { return *tree; }; + auto & Tree() { return *tree; }; }; } diff --git a/libsrc/gprim/geomobjects.hpp b/libsrc/gprim/geomobjects.hpp index ad215278..eb2f0030 100644 --- a/libsrc/gprim/geomobjects.hpp +++ b/libsrc/gprim/geomobjects.hpp @@ -270,7 +270,8 @@ namespace netgen for (int i = 0; i < D; i++) { if (p(i) < pmin(i)) pmin(i) = p(i); - else if (p(i) > pmax(i)) pmax(i) = p(i); + /* else */ if (p(i) > pmax(i)) pmax(i) = p(i); + // optimization invalid for empty-box ! } } diff --git a/libsrc/include/nginterface.h b/libsrc/include/nginterface.h index 3ddf24b8..24e90cc1 100644 --- a/libsrc/include/nginterface.h +++ b/libsrc/include/nginterface.h @@ -32,21 +32,25 @@ // max number of nodes per element -#define NG_ELEMENT_MAXPOINTS 12 +#define NG_ELEMENT_MAXPOINTS 20 // max number of nodes per surface element #define NG_SURFACE_ELEMENT_MAXPOINTS 8 +#ifndef PARALLEL + typedef int MPI_Comm; +#endif +namespace netgen { extern DLL_HEADER MPI_Comm ng_comm; } // implemented element types: enum NG_ELEMENT_TYPE { NG_PNT = 0, NG_SEGM = 1, NG_SEGM3 = 2, - NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, + NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, NG_QUAD8 = 14, NG_TET = 20, NG_TET10 = 21, - NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, - NG_HEX = 25 + NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, NG_PRISM15 = 27, NG_PYRAMID13 = 28, + NG_HEX = 25, NG_HEX20 = 26 }; typedef double NG_POINT[3]; // coordinates @@ -60,9 +64,9 @@ extern "C" { // load geometry from file DLL_HEADER void Ng_LoadGeometry (const char * filename); - + // load netgen mesh - DLL_HEADER void Ng_LoadMesh (const char * filename); + DLL_HEADER void Ng_LoadMesh (const char * filename, MPI_Comm comm = netgen::ng_comm); // load netgen mesh DLL_HEADER void Ng_LoadMeshFromString (const char * mesh_as_string); @@ -312,7 +316,7 @@ extern "C" { struct Ng_SolutionData { - string name; // name of gridfunction + std::string name; // name of gridfunction double * data; // solution values int components; // relevant (double) components in solution vector int dist; // # doubles per entry alignment! @@ -333,7 +337,7 @@ extern "C" { // redraw DLL_HEADER void Ng_Redraw(bool blocking = false); /// - DLL_HEADER void Ng_TclCmd(string cmd); + DLL_HEADER void Ng_TclCmd(std::string cmd); /// DLL_HEADER void Ng_SetMouseEventHandler (netgen::MouseEventHandler * handler); /// diff --git a/libsrc/include/nginterface_v2.hpp b/libsrc/include/nginterface_v2.hpp index 33a596db..bd432a43 100644 --- a/libsrc/include/nginterface_v2.hpp +++ b/libsrc/include/nginterface_v2.hpp @@ -12,9 +12,32 @@ C++ interface to Netgen */ +#ifndef NGINTERFACE + // implemented element types: +enum NG_ELEMENT_TYPE { + NG_PNT = 0, + NG_SEGM = 1, NG_SEGM3 = 2, + NG_TRIG = 10, NG_QUAD=11, NG_TRIG6 = 12, NG_QUAD6 = 13, NG_QUAD8 = 14, + NG_TET = 20, NG_TET10 = 21, + NG_PYRAMID = 22, NG_PRISM = 23, NG_PRISM12 = 24, NG_PRISM15 = 27, NG_PYRAMID13 = 28, + NG_HEX = 25, NG_HEX20 = 26 +}; + +enum NG_REFINEMENT_TYPE { NG_REFINE_H = 0, NG_REFINE_P = 1, NG_REFINE_HP = 2 }; +#endif + +#ifndef PARALLEL + typedef int MPI_Comm; +#endif + + namespace netgen { - + using namespace std; + using namespace ngcore; + + extern DLL_HEADER MPI_Comm ng_comm; + static constexpr int POINTINDEX_BASE = 1; struct T_EDGE2 @@ -44,6 +67,19 @@ namespace netgen size_t Size() const { return s; } T * Release() { T * hd = data; data = nullptr; return hd; } }; + + template + class Ng_BufferMS + { + size_t s; + T data[S]; + public: + Ng_BufferMS (size_t as) : s(as) { ; } + size_t Size() const { return s; } + T & operator[] (size_t i) { return data[i]; } + T operator[] (size_t i) const { return data[i]; } + }; + class Ng_Element { @@ -185,6 +221,7 @@ namespace netgen int operator[] (size_t i) const { return ptr[i]-POINTINDEX_BASE; } }; + /* class Ng_Edges { public: @@ -194,11 +231,11 @@ namespace netgen size_t Size() const { return ned; } int operator[] (size_t i) const { return ptr[i]-1; } }; - + */ public: Ng_Vertices vertices; - Ng_Edges edges; + // Ng_Edges edges; int surface_el; // -1 if face not on surface }; @@ -224,17 +261,24 @@ namespace netgen public: // Ngx_Mesh () { ; } // Ngx_Mesh(class Mesh * amesh) : mesh(amesh) { ; } - Ngx_Mesh(shared_ptr amesh = NULL); - void LoadMesh (const string & filename); - void LoadMesh (istream & str); + /** reuse a netgen-mesh **/ + Ngx_Mesh (shared_ptr amesh); + /** load a new mesh **/ + Ngx_Mesh (string filename, MPI_Comm acomm = netgen::ng_comm); + + void LoadMesh (const string & filename, MPI_Comm comm = netgen::ng_comm); + + void LoadMesh (istream & str, MPI_Comm comm = netgen::ng_comm); void SaveMesh (ostream & str) const; void UpdateTopology (); void DoArchive (Archive & archive); + MPI_Comm GetCommunicator() const; + virtual ~Ngx_Mesh(); - bool Valid () { return mesh != NULL; } + bool Valid () const { return mesh != NULL; } int GetDimension() const; int GetNLevels() const; @@ -282,7 +326,8 @@ namespace netgen template const Ng_Node GetNode (int nr) const; - + + Ng_BufferMS GetFaceEdges (int fnr) const; template int GetNNodes (); @@ -291,11 +336,17 @@ namespace netgen // 3D only // std::pair GetBoundaryNeighbouringDomains (int bnr); + template + void SetRefinementFlag (size_t elnr, bool flag); + void Curve (int order); + void Refine (NG_REFINEMENT_TYPE reftype, void (*taskmanager)(function) = &DummyTaskManager2, void (*tracer)(string, bool) = &DummyTracer2); + int GetHPElementLevel (int ei, int dir) const; + void GetParentNodes (int ni, int * parents) const; int GetParentElement (int ei) const; int GetParentSElement (int ei) const; @@ -312,13 +363,33 @@ namespace netgen int * const indices = NULL, int numind = 0) const; -#ifdef PARALLEL + // for MPI-parallel std::tuple GetDistantProcs (int nodetype, int locnum) const; -#endif shared_ptr GetMesh () const { return mesh; } shared_ptr SelectMesh () const; inline auto GetTimeStamp() const; + + + // also added from nginterface.h, still 1-based, need redesign + void HPRefinement (int levels, double parameter = 0.125, + bool setorders = true,bool ref_level = false); + size_t GetNP() const; + int GetSurfaceElementSurfaceNumber (size_t ei) const; + int GetSurfaceElementFDNumber (size_t ei) const; + + int GetElementOrder (int enr) const; + void GetElementOrders (int enr, int * ox, int * oy, int * oz) const; + void SetElementOrder (int enr, int order); + void SetElementOrders (int enr, int ox, int oy, int oz); + int GetSurfaceElementOrder (int enr) const; + void GetSurfaceElementOrders (int enr, int * ox, int * oy) const; + void SetSurfaceElementOrder (int enr, int order); + void SetSurfaceElementOrders (int enr, int ox, int oy); + int GetClusterRepVertex (int vi) const; + int GetClusterRepEdge (int edi) const; + int GetClusterRepFace (int fai) const; + int GetClusterRepElement (int eli) const; }; diff --git a/libsrc/include/nginterface_v2_impl.hpp b/libsrc/include/nginterface_v2_impl.hpp index 62e8068d..5b1efd6d 100644 --- a/libsrc/include/nginterface_v2_impl.hpp +++ b/libsrc/include/nginterface_v2_impl.hpp @@ -64,6 +64,13 @@ NGX_INLINE DLL_HEADER Ng_Element Ngx_Mesh :: GetElement<0> (size_t nr) const ret.facets.num = 1; ret.facets.base = 1; ret.facets.ptr = (int*)&el.pnum; + + if (mesh->GetDimension() == 1) + ret.mat = mesh->GetBCNamePtr(el.index-1); + else if (mesh->GetDimension() == 2) + ret.mat = mesh->GetCD2NamePtr(el.index-1); + else + ret.mat = mesh->GetCD3NamePtr(el.index-1); return ret; } diff --git a/libsrc/interface/CMakeLists.txt b/libsrc/interface/CMakeLists.txt index cb460fae..afd301d9 100644 --- a/libsrc/interface/CMakeLists.txt +++ b/libsrc/interface/CMakeLists.txt @@ -7,11 +7,9 @@ add_library(interface ${NG_LIB_TYPE} wuchemnitz.cpp writegmsh2.cpp writeOpenFOAM15x.cpp ) -if(NOT WIN32) - target_link_libraries(interface mesh csg geom2d) - target_link_libraries(interface visual) - install( TARGETS interface ${NG_INSTALL_DIR}) -endif(NOT WIN32) +target_link_libraries(interface mesh csg geom2d) +target_link_libraries(interface visual) +install( TARGETS interface ${NG_INSTALL_DIR}) install(FILES writeuser.hpp diff --git a/libsrc/interface/nginterface.cpp b/libsrc/interface/nginterface.cpp index ad0a72db..08be9458 100644 --- a/libsrc/interface/nginterface.cpp +++ b/libsrc/interface/nginterface.cpp @@ -117,14 +117,10 @@ void Ng_LoadMeshFromStream ( istream & input ) } - - -void Ng_LoadMesh (const char * filename) +void Ng_LoadMesh (const char * filename, MPI_Comm comm) { -#ifdef PARALLEL - MPI_Comm_size(MPI_COMM_WORLD, &ntasks); - MPI_Comm_rank(MPI_COMM_WORLD, &id); -#endif + int id = MyMPI_GetId(comm); + int ntasks = MyMPI_GetNTasks(comm); { ifstream infile(filename); @@ -134,11 +130,10 @@ void Ng_LoadMesh (const char * filename) if ( string(filename).find(".vol") == string::npos ) { -#ifdef PARALLEL if(ntasks>1) throw NgException("Not sure what to do with this?? Does this work with MPI??"); -#endif mesh.reset (new Mesh()); + mesh->SetCommunicator(comm); ReadFile(*mesh,filename); //mesh->SetGlobalH (mparam.maxh); //mesh->CalcLocalH(); @@ -146,12 +141,10 @@ void Ng_LoadMesh (const char * filename) } istream * infile; - char* buf; // for distributing geometry! + Array buf; // for distributing geometry! int strs; - #ifdef PARALLEL if( id == 0) { - #endif string fn(filename); if (fn.substr (fn.length()-3, 3) == ".gz") @@ -159,6 +152,7 @@ void Ng_LoadMesh (const char * filename) else infile = new ifstream (filename); mesh.reset (new Mesh()); + mesh->SetCommunicator(comm); mesh -> Load(*infile); SetGlobalMesh (mesh); @@ -168,12 +162,12 @@ void Ng_LoadMesh (const char * filename) geom_part << infile->rdbuf(); string geom_part_string = geom_part.str(); strs = geom_part_string.size(); - buf = new char[strs]; - memcpy(buf, geom_part_string.c_str(), strs*sizeof(char)); + // buf = new char[strs]; + buf.SetSize(strs); + memcpy(&buf[0], geom_part_string.c_str(), strs*sizeof(char)); } delete infile; -#ifdef PARALLEL if (ntasks > 1) { @@ -239,21 +233,21 @@ void Ng_LoadMesh (const char * filename) } // id==0 end else { mesh.reset (new Mesh()); + mesh->SetCommunicator(comm); SetGlobalMesh (mesh); mesh->SendRecvMesh(); } if(!ng_geometry && ntasks>1) { - /** Scatter the geometry-string **/ - MPI_Bcast(&strs, 1, MPI_INT, 0, MPI_COMM_WORLD); - if(id!=0) buf = new char[strs]; - MPI_Bcast(buf, strs, MPI_CHAR, 0, MPI_COMM_WORLD); - } +#ifdef PARALLEL + /** Scatter the geometry-string (no dummy-implementation in mpi_interface) **/ + MyMPI_Bcast(buf, comm); #endif + } if(!ng_geometry) { - infile = new istringstream(string((const char*)buf, (size_t)strs)); - delete[] buf; + infile = new istringstream(string((const char*)&buf[0], (size_t)strs)); + // delete[] buf; for (int i = 0; i < geometryregister.Size(); i++) { NetgenGeometry * hgeom = geometryregister[i]->LoadFromMeshFile (*infile); diff --git a/libsrc/interface/nginterface_v2.cpp b/libsrc/interface/nginterface_v2.cpp index fee63976..7fe34c1b 100644 --- a/libsrc/interface/nginterface_v2.cpp +++ b/libsrc/interface/nginterface_v2.cpp @@ -31,39 +31,39 @@ namespace netgen return hmesh; } - - Ngx_Mesh :: Ngx_Mesh (shared_ptr amesh) - { - if (amesh) - mesh = amesh; - else - mesh = netgen::mesh; - } + Ngx_Mesh :: Ngx_Mesh (shared_ptr amesh) + { mesh = amesh ? amesh : netgen::mesh; } + Ngx_Mesh :: Ngx_Mesh (string filename, MPI_Comm acomm) + { LoadMesh(filename, acomm); } - Ngx_Mesh * LoadMesh (const string & filename) + Ngx_Mesh * LoadMesh (const string & filename, MPI_Comm comm = netgen::ng_comm) { netgen::mesh.reset(); - Ng_LoadMesh (filename.c_str()); + Ng_LoadMesh (filename.c_str(), comm); return new Ngx_Mesh (netgen::mesh); } - void Ngx_Mesh :: LoadMesh (const string & filename) + void Ngx_Mesh :: LoadMesh (const string & filename, MPI_Comm comm) { netgen::mesh.reset(); - Ng_LoadMesh (filename.c_str()); + Ng_LoadMesh (filename.c_str(), comm); // mesh = move(netgen::mesh); mesh = netgen::mesh; } - void Ngx_Mesh :: LoadMesh (istream & ist) + void Ngx_Mesh :: LoadMesh (istream & ist, MPI_Comm comm) { netgen::mesh = make_shared(); + netgen::mesh->SetCommunicator(comm); netgen::mesh -> Load (ist); // mesh = move(netgen::mesh); mesh = netgen::mesh; SetGlobalMesh (mesh); } + MPI_Comm Ngx_Mesh :: GetCommunicator() const + { return Valid() ? mesh->GetCommunicator() : MPI_COMM_NULL; } + void Ngx_Mesh :: SaveMesh (ostream & ost) const { mesh -> Save (ost); @@ -71,7 +71,12 @@ namespace netgen void Ngx_Mesh :: DoArchive (Archive & archive) { - if (archive.Input()) mesh = make_shared(); +#ifdef PARALLEL + if (archive.Input()) { + mesh = make_shared(); + mesh->SetCommunicator(GetCommunicator()); + } +#endif mesh->DoArchive(archive); if (archive.Input()) { @@ -675,6 +680,37 @@ namespace netgen } + int Ngx_Mesh :: GetHPElementLevel (int ei, int dir) const + { + ei++; + int level = -1; + + if (mesh->hpelements) + { + int hpelnr = -1; + if (mesh->GetDimension() == 2) + hpelnr = mesh->SurfaceElement(ei).hp_elnr; + else + hpelnr = mesh->VolumeElement(ei).hp_elnr; + + if (hpelnr < 0) + throw NgException("Ngx_Mesh::GetHPElementLevel: Wrong hp-element number!"); + + if (dir == 1) + level = (*mesh->hpelements)[hpelnr].levelx; + else if (dir == 2) + level = (*mesh->hpelements)[hpelnr].levely; + else if (dir == 3) + level = (*mesh->hpelements)[hpelnr].levelz; + else + throw NgException("Ngx_Mesh::GetHPElementLevel: dir has to be 1, 2 or 3!"); + } + //else + // throw NgException("Ngx_Mesh::GetHPElementLevel only for HPRefinement implemented!"); + + return level; + } + int Ngx_Mesh :: GetParentElement (int ei) const { ei++; @@ -717,6 +753,16 @@ namespace netgen return mesh->GetIdentifications().GetType(idnr+1); } + Ng_BufferMS Ngx_Mesh::GetFaceEdges (int fnr) const + { + const MeshTopology & topology = mesh->GetTopology(); + ArrayMem ia; + topology.GetFaceEdges (fnr+1, ia); + Ng_BufferMS res(ia.Size()); + for (size_t i = 0; i < ia.Size(); i++) + res[i] = ia[i]-1; + return res; + } @@ -1042,6 +1088,19 @@ namespace netgen NgLock meshlock (mesh->MajorMutex(), true); mesh->BuildCurvedElements(order); } + + + template <> + DLL_HEADER void Ngx_Mesh :: SetRefinementFlag<2> (size_t elnr, bool flag) + { + mesh->SurfaceElement(elnr+1).SetRefinementFlag(flag); + } + + template <> + DLL_HEADER void Ngx_Mesh :: SetRefinementFlag<3> (size_t elnr, bool flag) + { + mesh->VolumeElement(elnr+1).SetRefinementFlag(flag); + } void Ngx_Mesh :: Refine (NG_REFINEMENT_TYPE reftype, void (*task_manager)(function), @@ -1070,13 +1129,127 @@ namespace netgen } + + + + // just copied with redesign + + size_t Ngx_Mesh::GetNP() const + { + return mesh->GetNP(); + } + + + int Ngx_Mesh::GetSurfaceElementSurfaceNumber (size_t ei) const + { + if (mesh->GetDimension() == 3) + return mesh->GetFaceDescriptor(mesh->SurfaceElement(ei).GetIndex()).SurfNr(); + else + return mesh->LineSegment(ei).si; + } + int Ngx_Mesh::GetSurfaceElementFDNumber (size_t ei) const + { + if (mesh->GetDimension() == 3) + return mesh->SurfaceElement(ei).GetIndex(); + else + return -1; + } + + + void Ngx_Mesh::HPRefinement (int levels, double parameter, bool setorders, + bool ref_level) + { + NgLock meshlock (mesh->MajorMutex(), true); + Refinement & ref = const_cast (mesh->GetGeometry()->GetRefinement()); + ::netgen::HPRefinement (*mesh, &ref, levels, parameter, setorders, ref_level); + } + +int Ngx_Mesh::GetElementOrder (int enr) const +{ + if (mesh->GetDimension() == 3) + return mesh->VolumeElement(enr).GetOrder(); + else + return mesh->SurfaceElement(enr).GetOrder(); +} + +void Ngx_Mesh::GetElementOrders (int enr, int * ox, int * oy, int * oz) const +{ + if (mesh->GetDimension() == 3) + mesh->VolumeElement(enr).GetOrder(*ox, *oy, *oz); + else + mesh->SurfaceElement(enr).GetOrder(*ox, *oy, *oz); +} + +void Ngx_Mesh::SetElementOrder (int enr, int order) +{ + if (mesh->GetDimension() == 3) + return mesh->VolumeElement(enr).SetOrder(order); + else + return mesh->SurfaceElement(enr).SetOrder(order); +} + +void Ngx_Mesh::SetElementOrders (int enr, int ox, int oy, int oz) +{ + if (mesh->GetDimension() == 3) + mesh->VolumeElement(enr).SetOrder(ox, oy, oz); + else + mesh->SurfaceElement(enr).SetOrder(ox, oy); +} + + +int Ngx_Mesh::GetSurfaceElementOrder (int enr) const +{ + return mesh->SurfaceElement(enr).GetOrder(); +} + +int Ngx_Mesh::GetClusterRepVertex (int pi) const +{ + return mesh->GetClusters().GetVertexRepresentant(pi); +} + +int Ngx_Mesh::GetClusterRepEdge (int pi) const +{ + return mesh->GetClusters().GetEdgeRepresentant(pi); +} + +int Ngx_Mesh::GetClusterRepFace (int pi) const +{ + return mesh->GetClusters().GetFaceRepresentant(pi); +} + +int Ngx_Mesh::GetClusterRepElement (int pi) const +{ + return mesh->GetClusters().GetElementRepresentant(pi); +} + + + + +//HERBERT: falsche Anzahl von Argumenten +//void Ngx_Mesh::GetSurfaceElementOrders (int enr, int * ox, int * oy, int * oz) +void Ngx_Mesh::GetSurfaceElementOrders (int enr, int * ox, int * oy) const +{ + int d; + mesh->SurfaceElement(enr).GetOrder(*ox, *oy, d); +} + +void Ngx_Mesh::SetSurfaceElementOrder (int enr, int order) +{ + return mesh->SurfaceElement(enr).SetOrder(order); +} + +void Ngx_Mesh::SetSurfaceElementOrders (int enr, int ox, int oy) +{ + mesh->SurfaceElement(enr).SetOrder(ox, oy); +} + -#ifdef PARALLEL + std::tuple Ngx_Mesh :: GetDistantProcs (int nodetype, int locnum) const { - +#ifdef PARALLEL switch (nodetype) { case 0: @@ -1097,10 +1270,10 @@ namespace netgen default: return std::tuple(0,nullptr); } - } - +#else + return std::tuple(0,nullptr); #endif - + } } diff --git a/libsrc/linalg/CMakeLists.txt b/libsrc/linalg/CMakeLists.txt index a5a5a4f2..a8ab1b5a 100644 --- a/libsrc/linalg/CMakeLists.txt +++ b/libsrc/linalg/CMakeLists.txt @@ -1,8 +1,8 @@ -add_library( la OBJECT - densemat.cpp polynomial.cpp bfgs.cpp linopt.cpp linsearch.cpp - ) - -set_target_properties(la PROPERTIES POSITION_INDEPENDENT_CODE ON ) +add_library( la INTERFACE ) +set(sdir ${CMAKE_CURRENT_SOURCE_DIR}) +target_sources( la INTERFACE + ${sdir}/densemat.cpp ${sdir}/polynomial.cpp ${sdir}/bfgs.cpp ${sdir}/linopt.cpp ${sdir}/linsearch.cpp +) install(FILES densemat.hpp linalg.hpp opti.hpp diff --git a/libsrc/meshing/CMakeLists.txt b/libsrc/meshing/CMakeLists.txt index bc9f7f18..12fe70ae 100644 --- a/libsrc/meshing/CMakeLists.txt +++ b/libsrc/meshing/CMakeLists.txt @@ -1,12 +1,4 @@ add_definitions(-DNGINTERFACE_EXPORTS) -if(NOT WIN32) - set(mesh_object_libs - $ - $ - $ - ) -endif(NOT WIN32) - add_library(mesh ${NG_LIB_TYPE} adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp clusters.cpp curvedelems.cpp delaunay.cpp delaunay2d.cpp @@ -29,10 +21,10 @@ if(APPLE) set_target_properties( mesh PROPERTIES SUFFIX ".so") endif(APPLE) -if(NOT WIN32) - target_link_libraries( mesh ngcore ${ZLIB_LIBRARIES} ${MPI_CXX_LIBRARIES} ${PYTHON_LIBRARIES} ${METIS_LIBRARY}) - install( TARGETS mesh ${NG_INSTALL_DIR}) -endif(NOT WIN32) +target_link_libraries( mesh PUBLIC ngcore PRIVATE gprim la gen ) + +target_link_libraries( mesh PUBLIC ${ZLIB_LIBRARIES} ${MPI_CXX_LIBRARIES} ${PYTHON_LIBRARIES} ${METIS_LIBRARY}) +install( TARGETS mesh ${NG_INSTALL_DIR}) install(FILES adfront2.hpp adfront3.hpp basegeom.hpp bcfunctions.hpp bisect.hpp diff --git a/libsrc/meshing/adfront3.hpp b/libsrc/meshing/adfront3.hpp index 738de122..9b7f3818 100644 --- a/libsrc/meshing/adfront3.hpp +++ b/libsrc/meshing/adfront3.hpp @@ -233,6 +233,7 @@ public: /// const MiniElement2d & GetFace (int i) const { return faces.Get(i).Face(); } + const auto & Faces() const { return faces; } /// void Print () const; /// diff --git a/libsrc/meshing/bisect.cpp b/libsrc/meshing/bisect.cpp index 9edf90d2..d73718be 100644 --- a/libsrc/meshing/bisect.cpp +++ b/libsrc/meshing/bisect.cpp @@ -3398,6 +3398,7 @@ namespace netgen size_t nsel = mtris.Size(); NgProfiler::StartTimer (timer_bisecttrig); + (*opt.tracer)("Bisect trigs", false); for (size_t i = 0; i < nsel; i++) if (mtris[i].marked) { @@ -3440,7 +3441,7 @@ namespace netgen } NgProfiler::StopTimer (timer_bisecttrig); - + (*opt.tracer)("Bisect trigs", true); int nquad = mquads.Size(); for (int i = 1; i <= nquad; i++) diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp index 96a944fc..182945d3 100644 --- a/libsrc/meshing/curvedelems.cpp +++ b/libsrc/meshing/curvedelems.cpp @@ -553,14 +553,18 @@ namespace netgen order = 1; + MPI_Comm curve_comm; #ifdef PARALLEL enum { MPI_TAG_CURVE = MPI_TAG_MESH+20 }; const ParallelMeshTopology & partop = mesh.GetParallelTopology (); - MPI_Comm curve_comm; - MPI_Comm_dup (MPI_COMM_WORLD, &curve_comm); + MPI_Comm_dup (mesh.GetCommunicator(), &curve_comm); Array procs; +#else + curve_comm = ng_comm; // dummy! #endif + int rank = MyMPI_GetId(curve_comm); + int ntasks = MyMPI_GetNTasks(curve_comm); if (working) order = aorder; @@ -1708,6 +1712,7 @@ namespace netgen case TRIG : info.nv = 3; break; case QUAD : info.nv = 4; break; case TRIG6: info.nv = 6; break; + case QUAD8 : info.nv = 8; break; default: cerr << "undef element in CalcSurfaceTrafo" << endl; } @@ -1927,6 +1932,24 @@ namespace netgen break; } + + case QUAD8: + { + auto x = xi(0), y = xi(1); + shapes(0) = (1-x)*(1-y); + shapes(1) = x*(1-y); + shapes(2) = x*y; + shapes(3) = (1-x)*y; + shapes(4) = 4*(1-x)*x*(1-y); + shapes(5) = 4*(1-x)*x*y; + shapes(6) = 4*(1-y)*y*(1-x); + shapes(7) = 4*(1-y)*y*x; + shapes(0) -= 0.5*(shapes(4)+shapes(6)); + shapes(1) -= 0.5*(shapes(4)+shapes(7)); + shapes(2) -= 0.5*(shapes(5)+shapes(7)); + shapes(3) -= 0.5*(shapes(5)+shapes(6)); + break; + } default: throw NgException("CurvedElements::CalcShape 2d, element type not handled"); @@ -2266,7 +2289,7 @@ namespace netgen } break; } - case QUAD: + case QUAD: { if (info.order >= 2) return false; // not yet supported AutoDiff<2,T> lami[4] = { (1-x)*(1-y), x*(1-y), x*y, (1-x)*y }; @@ -2278,6 +2301,33 @@ namespace netgen } break; } + case QUAD8: + { + // AutoDiff<2,T> lami[4] = { (1-x)*(1-y), x*(1-y), x*y, (1-x)*y }; + AutoDiff<2,T> lami[8] = + { (1-x)*(1-y), + x*(1-y), + x*y, + (1-x)*y, + 4*(1-x)*x*(1-y), + 4*(1-x)*x*y, + 4*(1-y)*y*(1-x), + 4*(1-y)*y*x }; + + lami[0] -= 0.5*(lami[4]+lami[6]); + lami[1] -= 0.5*(lami[4]+lami[7]); + lami[2] -= 0.5*(lami[5]+lami[7]); + lami[3] -= 0.5*(lami[5]+lami[6]); + + for (int j = 0; j < 8; j++) + { + Point<3> p = mesh[el[j]]; + for (int k = 0; k < DIM_SPACE; k++) + mapped_x[k] += p(k) * lami[j]; + } + break; + } + default: return false; } @@ -2740,6 +2790,32 @@ namespace netgen break; } + case PRISM15: + { + shapes = 0.0; + T x = xi(0); + T y = xi(1); + T z = xi(2); + T lam = 1-x-y; + T lamz = 1-z; + shapes[0] = (2*x*x-x) * (2*lamz*lamz-lamz); + shapes[1] = (2*y*y-y) * (2*lamz*lamz-lamz); + shapes[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz); + shapes[3] = (2*x*x-x) * (2*z*z-z); + shapes[4] = (2*y*y-y) * (2*z*z-z); + shapes[5] = (2*lam*lam-lam) * (2*z*z-z); + shapes[6] = 4 * x * y * (2*lamz*lamz-lamz); + shapes[7] = 4 * x * lam * (2*lamz*lamz-lamz); + shapes[8] = 4 * y * lam * (2*lamz*lamz-lamz); + shapes[9] = x * 4 * z * (1-z); + shapes[10] = y * 4 * z * (1-z); + shapes[11] = lam * 4 * z * (1-z); + shapes[12] = 4 * x * y * (2*z*z-z); + shapes[13] = 4 * x * lam * (2*z*z-z); + shapes[14] = 4 * y * lam * (2*z*z-z); + break; + } + case PYRAMID: { shapes = 0.0; @@ -2789,6 +2865,29 @@ namespace netgen break; } + case PYRAMID13: + { + shapes = 0.0; + T x = xi(0); + T y = xi(1); + T z = xi(2); + z *= 1-1e-12; + shapes[0] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (-2*x - z + 2)*(-2*y - z + 2))*(-0.5*x - 0.5*y - 0.5*z + 0.25); + shapes[1] = (0.5*x - 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(-2*y - z + 2)); + shapes[2] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(2*y + z))*(0.5*x + 0.5*y + 0.5*z - 0.75); + shapes[3] = (-0.5*x + 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*y + z)*(-2*x - z + 2)); + shapes[4] = z*(2*z - 1); + shapes[5] = 2*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shapes[6] = 4*x*y*(-2*x - 2*z + 2)/(-2*z + 2); + shapes[7] = 2*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shapes[8] = 4*x*y*(-2*y - 2*z + 2)/(-2*z + 2); + shapes[9] = z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1); + shapes[10] = 2*x*z*(-2*y - 2*z + 2)/(-z + 1); + shapes[11] = 4*x*y*z/(-z + 1); + shapes[12] = 2*y*z*(-2*x - 2*z + 2)/(-z + 1); + break; + } + case HEX: { shapes = 0.0; @@ -2839,6 +2938,46 @@ namespace netgen break; + } + + case HEX20: + { + shapes = 0.0; + T x = xi(0); + T y = xi(1); + T z = xi(2); + + shapes[0] = (1-x)*(1-y)*(1-z); + shapes[1] = x *(1-y)*(1-z); + shapes[2] = x * y *(1-z); + shapes[3] = (1-x)* y *(1-z); + shapes[4] = (1-x)*(1-y)*(z); + shapes[5] = x *(1-y)*(z); + shapes[6] = x * y *(z); + shapes[7] = (1-x)* y *(z); + + T sigma[8]={(1-x)+(1-y)+(1-z),x+(1-y)+(1-z),x+y+(1-z),(1-x)+y+(1-z), + (1-x)+(1-y)+z,x+(1-y)+z,x+y+z,(1-x)+y+z}; + + static const int e[12][2] = + { + { 0, 1 }, { 2, 3 }, { 3, 0 }, { 1, 2 }, + { 4, 5 }, { 6, 7 }, { 7, 4 }, { 5, 6 }, + { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }, + }; + + for (int i = 0; i < 12; i++) + { + T lame = shapes[e[i][0]]+shapes[e[i][1]]; + T xi = sigma[e[i][1]]-sigma[e[i][0]]; + shapes[8+i] = (1-xi*xi)*lame; + } + for (int i = 0; i < 12; i++) + { + shapes[e[i][0]] -= 0.5 * shapes[8+i]; + shapes[e[i][1]] -= 0.5 * shapes[8+i]; + } + break; } default: @@ -3210,6 +3349,36 @@ namespace netgen } + case PRISM15: + { + AutoDiff<3,T> x(xi(0), 0); + AutoDiff<3,T> y(xi(1), 1); + AutoDiff<3,T> z(xi(2), 2); + AutoDiff<3,T> ad[15]; + AutoDiff<3,T> lam = 1-x-y; + AutoDiff<3,T> lamz = 1-z; + + ad[0] = (2*x*x-x) * (2*lamz*lamz-lamz); + ad[1] = (2*y*y-y) * (2*lamz*lamz-lamz); + ad[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz); + ad[3] = (2*x*x-x) * (2*z*z-z); + ad[4] = (2*y*y-y) * (2*z*z-z); + ad[5] = (2*lam*lam-lam) * (2*z*z-z); + ad[6] = 4 * x * y * (2*lamz*lamz-lamz); + ad[7] = 4 * x * lam * (2*lamz*lamz-lamz); + ad[8] = 4 * y * lam * (2*lamz*lamz-lamz); + ad[9] = x * 4 * z * (1-z); + ad[10] = y * 4 * z * (1-z); + ad[11] = lam * 4 * z * (1-z); + ad[12] = 4 * x * y * (2*z*z-z); + ad[13] = 4 * x * lam * (2*z*z-z); + ad[14] = 4 * y * lam * (2*z*z-z); + + for(int i=0; i<15; i++) + for(int j=0; j<3; j++) + dshapes(i,j) = ad[i].DValue(j); + break; + } case PYRAMID: { // if (typeid(T) == typeid(SIMD)) return; @@ -3307,6 +3476,54 @@ namespace netgen break; } + case PYRAMID13: + { + T x = xi(0); + T y = xi(1); + T z = xi(2); + z *= 1-1e-12; + dshapes(0,0) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(4*y + 2*z + 2*z*(2*y + z - 1)/(-z + 1) - 4); + dshapes(0,1) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(4*x + 2*z + 2*z*(2*x + z - 1)/(-z + 1) - 4); + dshapes(0,2) = 0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(-2*x - z + 2)*(-2*y - z + 2) + (-0.5*x - 0.5*y - 0.5*z + 0.25)*(2*x + 2*y + 2*z + z*(2*x + z - 1)/(-z + 1) + z*(2*y + z - 1)/(-z + 1) + z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) - 5 + (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(1,0) = -0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(-2*y - z + 2) + (0.5*x - 0.5*y - 0.25)*(-4*y - 2*z - 2*z*(2*y + z - 1)/(-z + 1) + 4); + dshapes(1,1) = 0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(2*x + z)*(-2*y - z + 2) + (-4*x - 2*z - 2*z*(2*x + z - 1)/(-z + 1))*(0.5*x - 0.5*y - 0.25); + dshapes(1,2) = (0.5*x - 0.5*y - 0.25)*(-2*x - 2*y - 2*z - z*(2*x + z - 1)/(-z + 1) - z*(2*y + z - 1)/(-z + 1) - z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) + 1 - (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(2,0) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (4*y + 2*z + 2*z*(2*y + z - 1)/(-z + 1))*(0.5*x + 0.5*y + 0.5*z - 0.75); + dshapes(2,1) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (4*x + 2*z + 2*z*(2*x + z - 1)/(-z + 1))*(0.5*x + 0.5*y + 0.5*z - 0.75); + dshapes(2,2) = -0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*x + z)*(2*y + z) + (0.5*x + 0.5*y + 0.5*z - 0.75)*(2*x + 2*y + 2*z + z*(2*x + z - 1)/(-z + 1) + z*(2*y + z - 1)/(-z + 1) + z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) - 1 + (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(3,0) = 0.5*z + 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) - 0.5*(2*y + z)*(-2*x - z + 2) + (-0.5*x + 0.5*y - 0.25)*(-4*y - 2*z - 2*z*(2*y + z - 1)/(-z + 1)); + dshapes(3,1) = -0.5*z - 0.5*z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + 0.5*(2*y + z)*(-2*x - z + 2) + (-0.5*x + 0.5*y - 0.25)*(-4*x - 2*z - 2*z*(2*x + z - 1)/(-z + 1) + 4); + dshapes(3,2) = (-0.5*x + 0.5*y - 0.25)*(-2*x - 2*y - 2*z - z*(2*x + z - 1)/(-z + 1) - z*(2*y + z - 1)/(-z + 1) - z*(2*x + z - 1)*(2*y + z - 1)/((-z + 1)*(-z + 1)) + 1 - (2*x + z - 1)*(2*y + z - 1)/(-z + 1)); + dshapes(4,0) = 0; + dshapes(4,1) = 0; + dshapes(4,2) = 4*z - 1; + dshapes(5,0) = -4*x*(-2*y - 2*z + 2)/(-2*z + 2) + 2*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(5,1) = -4*x*(-2*x - 2*z + 2)/(-2*z + 2); + dshapes(5,2) = -4*x*(-2*x - 2*z + 2)/(-2*z + 2) - 4*x*(-2*y - 2*z + 2)/(-2*z + 2) + 4*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(6,0) = -8*x*y/(-2*z + 2) + 4*y*(-2*x - 2*z + 2)/(-2*z + 2); + dshapes(6,1) = 4*x*(-2*x - 2*z + 2)/(-2*z + 2); + dshapes(6,2) = -8*x*y/(-2*z + 2) + 8*x*y*(-2*x - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(7,0) = -4*y*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(7,1) = -4*y*(-2*x - 2*z + 2)/(-2*z + 2) + 2*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(7,2) = -4*y*(-2*x - 2*z + 2)/(-2*z + 2) - 4*y*(-2*y - 2*z + 2)/(-2*z + 2) + 4*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(8,0) = 4*y*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(8,1) = -8*x*y/(-2*z + 2) + 4*x*(-2*y - 2*z + 2)/(-2*z + 2); + dshapes(8,2) = -8*x*y/(-2*z + 2) + 8*x*y*(-2*y - 2*z + 2)/((-2*z + 2)*(-2*z + 2)); + dshapes(9,0) = -2*z*(-2*y - 2*z + 2)/(-z + 1); + dshapes(9,1) = -2*z*(-2*x - 2*z + 2)/(-z + 1); + dshapes(9,2) = -2*z*(-2*x - 2*z + 2)/(-z + 1) - 2*z*(-2*y - 2*z + 2)/(-z + 1) + z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/((-z + 1)*(-z + 1)) + (-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1); + dshapes(10,0) = 2*z*(-2*y - 2*z + 2)/(-z + 1); + dshapes(10,1) = -4*x*z/(-z + 1); + dshapes(10,2) = -4*x*z/(-z + 1) + 2*x*z*(-2*y - 2*z + 2)/((-z + 1)*(-z + 1)) + 2*x*(-2*y - 2*z + 2)/(-z + 1); + dshapes(11,0) = 4*y*z/(-z + 1); + dshapes(11,1) = 4*x*z/(-z + 1); + dshapes(11,2) = 4*x*y*z/((-z + 1)*(-z + 1)) + 4*x*y/(-z + 1); + dshapes(12,0) = -4*y*z/(-z + 1); + dshapes(12,1) = 2*z*(-2*x - 2*z + 2)/(-z + 1); + dshapes(12,2) = -4*y*z/(-z + 1) + 2*y*z*(-2*x - 2*z + 2)/((-z + 1)*(-z + 1)) + 2*y*(-2*x - 2*z + 2)/(-z + 1); + break; + } + case HEX: { // if (typeid(T) == typeid(SIMD)) return; @@ -3443,12 +3660,49 @@ namespace netgen *testout << "quad, num dshape = " << endl << dshapes << endl; */ break; - - - - break; } - + case HEX20: + { + AutoDiff<3,T> x(xi(0), 0); + AutoDiff<3,T> y(xi(1), 1); + AutoDiff<3,T> z(xi(2), 2); + AutoDiff<3,T> ad[20]; + + ad[0] = (1-x)*(1-y)*(1-z); + ad[1] = x *(1-y)*(1-z); + ad[2] = x * y *(1-z); + ad[3] = (1-x)* y *(1-z); + ad[4] = (1-x)*(1-y)*(z); + ad[5] = x *(1-y)*(z); + ad[6] = x * y *(z); + ad[7] = (1-x)* y *(z); + + AutoDiff<3,T> sigma[8]={(1-x)+(1-y)+(1-z),x+(1-y)+(1-z),x+y+(1-z),(1-x)+y+(1-z), + (1-x)+(1-y)+z,x+(1-y)+z,x+y+z,(1-x)+y+z}; + + static const int e[12][2] = + { + { 0, 1 }, { 2, 3 }, { 3, 0 }, { 1, 2 }, + { 4, 5 }, { 6, 7 }, { 7, 4 }, { 5, 6 }, + { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }, + }; + + for (int i = 0; i < 12; i++) + { + auto lame = ad[e[i][0]]+ad[e[i][1]]; + auto xi = sigma[e[i][1]]-sigma[e[i][0]]; + ad[8+i] = (1-xi*xi)*lame; + } + for (int i = 0; i < 12; i++) + { + ad[e[i][0]] -= 0.5 * ad[8+i]; + ad[e[i][1]] -= 0.5 * ad[8+i]; + } + for (int i = 0; i < 20; i++) + for (int j = 0; j < 3; j++) + dshapes(i,j) = ad[i].DValue(j); + break; + } default: throw NgException("CurvedElements::CalcDShape 3d, element type not handled"); } @@ -3854,8 +4108,9 @@ namespace netgen case TRIG : info.nv = 3; break; case QUAD : info.nv = 4; break; case TRIG6: info.nv = 6; break; + case QUAD8 : info.nv = 8; break; default: - cerr << "undef element in CalcMultPointSurfaceTrao" << endl; + cerr << "undef element in CalcMultPointSurfaceTrafo" << endl; } info.ndof = info.nv; diff --git a/libsrc/meshing/delaunay.cpp b/libsrc/meshing/delaunay.cpp index bf0a5174..cbd2a73b 100644 --- a/libsrc/meshing/delaunay.cpp +++ b/libsrc/meshing/delaunay.cpp @@ -6,7 +6,6 @@ namespace netgen { - static const int deltetfaces[][3] = { { 1, 2, 3 }, { 2, 0, 3 }, @@ -14,10 +13,6 @@ namespace netgen { 1, 0, 2 } }; - - - - class DelaunayTet { PointIndex pnums[4]; @@ -41,9 +36,6 @@ namespace netgen PointIndex & operator[] (int i) { return pnums[i]; } PointIndex operator[] (int i) const { return pnums[i]; } - int & NB1(int i) { return nb[i-1]; } - int NB1(int i) const { return nb[i-1]; } - int & NB(int i) { return nb[i]; } int NB(int i) const { return nb[i]; } @@ -57,27 +49,6 @@ namespace netgen return i; return 3; } - - void GetFace1 (int i, INDEX_3 & face) const - { - face.I(1) = pnums[deltetfaces[i-1][0]]; - face.I(2) = pnums[deltetfaces[i-1][1]]; - face.I(3) = pnums[deltetfaces[i-1][2]]; - } - - void GetFace (int i, INDEX_3 & face) const - { - face.I(1) = pnums[deltetfaces[i][0]]; - face.I(2) = pnums[deltetfaces[i][1]]; - face.I(3) = pnums[deltetfaces[i][2]]; - } - - INDEX_3 GetFace1 (int i) const - { - return INDEX_3 (pnums[deltetfaces[i-1][0]], - pnums[deltetfaces[i-1][1]], - pnums[deltetfaces[i-1][2]]); - } INDEX_3 GetFace (int i) const { @@ -85,13 +56,13 @@ namespace netgen pnums[deltetfaces[i][1]], pnums[deltetfaces[i][2]]); } - - void GetFace1 (int i, Element2d & face) const + + void GetFace (int i, Element2d & face) const { // face.SetType(TRIG); - face[0] = pnums[deltetfaces[i-1][0]]; - face[1] = pnums[deltetfaces[i-1][1]]; - face[2] = pnums[deltetfaces[i-1][2]]; + face[0] = pnums[deltetfaces[i][0]]; + face[1] = pnums[deltetfaces[i][1]]; + face[2] = pnums[deltetfaces[i][2]]; } }; @@ -101,8 +72,6 @@ namespace netgen - - /* Table to maintain neighbour elements */ @@ -135,7 +104,7 @@ namespace netgen // get neighbour of element elnr in direction fnr int GetNB (int elnr, int fnr) { - return tets.Get(elnr).NB1(fnr); + return tets.Get(elnr).NB(fnr); } // @@ -181,7 +150,6 @@ namespace netgen - /* connected lists of cosphereical elements */ @@ -253,24 +221,23 @@ namespace netgen Array & freelist, SphereList & list, IndexSet & insphere, IndexSet & closesphere) { + // static Timer t("Meshing3::AddDelaunayPoint"); RegionTimer reg(t); /* find any sphere, such that newp is contained in */ - DelaunayTet el; int cfelind = -1; const Point<3> * pp[4]; Point<3> pc; - double r2; Point3d tpmin, tpmax; tettree.GetIntersecting (newp, newp, treesearch); + double quot,minquot(1e20); - for (int j = 0; j < treesearch.Size(); j++) + for (auto jjj : treesearch) { - int jjj = treesearch[j]; quot = Dist2 (centers.Get(jjj), newp) / radi2.Get(jjj); if((cfelind == -1 || quot < 0.99*minquot) && quot < 1) @@ -284,46 +251,12 @@ namespace netgen } - /* - int i, j, k, l; - if (!felind) - { - cerr << "not in any sphere, 1" << endl; - // old, non tree search - - double mindist = 1e10; - for (j = 1; j <= tempels.Size(); j++) - { - if (tempels.Get(j).PNum(1)) - { - double toofar = - Dist2 (centers.Get(j), newp) - radi2.Get(j); - if (toofar < mindist || toofar < 1e-7) - { - mindist = toofar; - cout << " dist2 = " << Dist2 (centers.Get(j), newp) - << " radi2 = " << radi2.Get(j) << endl; - } - if (toofar < 0) - { - el = tempels.Get(j); - felind = j; - cout << "sphere found !" << endl; - break; - } - } - } - cout << "point is too far from sheres: " << mindist << endl; - } - */ - if (cfelind == -1) { PrintWarning ("Delaunay, point not in any sphere"); return; } - /* insphere: point is in sphere -> delete element closesphere: point is close to sphere -> considered for same center @@ -376,7 +309,7 @@ namespace netgen // check neighbour-tets for (int j = starti; j < nstarti; j++) - for (int k = 1; k <= 4; k++) + for (int k = 0; k < 4; k++) { int helind = insphere.GetArray().Get(j); int nbind = meshnb.GetNB (helind, k); @@ -404,61 +337,48 @@ namespace netgen } else { - /* - Element2d face; - tempels.Get(helind).GetFace (k, face); + INDEX_3 i3 = tempels.Get(helind).GetFace (k); - const Point3d & p1 = mesh.Point (face.PNum(1)); - const Point3d & p2 = mesh.Point (face[1]); - const Point3d & p3 = mesh.Point (face[2]); - */ + const Point<3> & p1 = mesh.Point ( PointIndex (i3.I1()) ); + const Point<3> & p2 = mesh.Point ( PointIndex (i3.I2()) ); + const Point<3> & p3 = mesh.Point ( PointIndex (i3.I3()) ); - INDEX_3 i3 = tempels.Get(helind).GetFace (k-1); - - const Point3d & p1 = mesh.Point ( PointIndex (i3.I1()) ); - const Point3d & p2 = mesh.Point ( PointIndex (i3.I2()) ); - const Point3d & p3 = mesh.Point ( PointIndex (i3.I3()) ); - - - Vec3d v1(p1, p2); - Vec3d v2(p1, p3); - Vec3d n = Cross (v1, v2); + Vec<3> v1 = p2-p1; + Vec<3> v2 = p3-p1; + Vec<3> n = Cross (v1, v2); n /= n.Length(); - if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k-1])) > 0) + if (n * Vec3d (p1, mesh.Point (tempels.Get(helind)[k])) > 0) n *= -1; double dist = n * Vec3d (p1, newp); - if (dist > -1e-10) // 1e-10 { insphere.Add (nbind); changed = 1; } - - } } } } // while (changed) - // (*testout) << "newels: " << endl; - Array newels; + + // Array newels; + Array newels; Element2d face(TRIG); - for (int j = 1; j <= insphere.GetArray().Size(); j++) - for (int k = 1; k <= 4; k++) + for (int celind : insphere.GetArray()) + for (int k = 0; k < 4; k++) { - // int elind = insphere.GetArray().Get(j); - int celind = insphere.GetArray().Get(j); int nbind = meshnb.GetNB (celind, k); if (!nbind || !insphere.IsIn (nbind)) { - tempels.Get (celind).GetFace1 (k, face); + tempels.Get (celind).GetFace (k, face); - Element newel(TET); + // Element newel(TET); + DelaunayTet newel; for (int l = 0; l < 3; l++) newel[l] = face[l]; newel[3] = newpi; @@ -471,7 +391,7 @@ namespace netgen n.Normalize(); if (n * Vec3d(mesh.Point (face[0]), - mesh.Point (tempels.Get(insphere.GetArray().Get(j))[k-1])) + mesh.Point (tempels.Get(celind)[k])) > 0) n *= -1; @@ -492,38 +412,36 @@ namespace netgen meshnb.ResetFaceHT (10*insphere.GetArray().Size()+1); - for (int j = 1; j <= insphere.GetArray().Size(); j++) + for (auto celind : insphere.GetArray()) { - // int elind = - int celind = insphere.GetArray().Get(j); - meshnb.Delete (celind); list.DeleteElement (celind); for (int k = 0; k < 4; k++) tempels.Elem(celind)[k] = -1; - // ((ADTree6&)tettree.Tree()).DeleteElement (celind); tettree.DeleteElement (celind); freelist.Append (celind); } - int hasclose = 0; - for (int j = 1; j <= closesphere.GetArray().Size(); j++) + + + bool hasclose = false; + for (int ind : closesphere.GetArray()) { - int ind = closesphere.GetArray().Get(j); if (!insphere.IsIn(ind) && fabs (Dist2 (centers.Get (ind), newp) - radi2.Get(ind)) < 1e-8 ) - hasclose = 1; + hasclose = true; } for (int j = 1; j <= newels.Size(); j++) { + const auto & newel = newels.Get(j); int nelind; if (!freelist.Size()) { - tempels.Append (newels.Get(j)); + tempels.Append (newel); nelind = tempels.Size(); } else @@ -531,29 +449,29 @@ namespace netgen nelind = freelist.Last(); freelist.DeleteLast(); - tempels.Elem(nelind) = newels.Get(j); + tempels.Elem(nelind) = newel; } meshnb.Add (nelind); list.AddElement (nelind); for (int k = 0; k < 4; k++) - pp[k] = &mesh.Point (newels.Get(j)[k]); + pp[k] = &mesh.Point (newel[k]); if (CalcSphereCenter (&pp[0], pc) ) { PrintSysError ("Delaunay: New tet is flat"); (*testout) << "new tet is flat" << endl; - for (int k = 1; k <= 4; k++) - (*testout) << newels.Get(j).PNum(k) << " "; + for (int k = 0; k < 4; k++) + (*testout) << newel[k] << " "; (*testout) << endl; - for (int k = 1; k <= 4; k++) + for (int k = 0; k < 4; k++) (*testout) << *pp[k-1] << " "; (*testout) << endl; } - r2 = Dist2 (*pp[0], pc); + double r2 = Dist2 (*pp[0], pc); if (hasclose) for (int k = 1; k <= closesphere.GetArray().Size(); k++) { @@ -602,50 +520,41 @@ namespace netgen Array & tempels, int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax) { - Array > centers; + static Timer t("Meshing3::Delaunay1"); RegionTimer reg(t); + static Timer tloop("Meshing3::Delaunay1 loop"); + + Array> centers; Array radi2; - Point3d tpmin, tpmax; + Box<3> bbox(Box<3>::EMPTY_BOX); + for (auto & face : adfront->Faces()) + for (PointIndex pi : face.Face().PNums()) + bbox.Add (mesh.Point(pi)); - // new: local box - mesh.GetBox (pmax, pmin); // lower bound for pmax, upper for pmin - for (int i = 1; i <= adfront->GetNF(); i++) - { - const MiniElement2d & face = adfront->GetFace(i); - for (PointIndex pi : face.PNums()) - { - pmin.SetToMin (mesh.Point (pi)); - pmax.SetToMax (mesh.Point (pi)); - } - } - for (PointIndex pi : mesh.LockedPoints()) - { - pmin.SetToMin (mesh.Point (pi)); - pmax.SetToMax (mesh.Point (pi)); - } - + bbox.Add (mesh.Point (pi)); + pmin = bbox.PMin(); + pmax = bbox.PMax(); - Vec3d vdiag(pmin, pmax); + + Vec<3> vdiag = pmax-pmin; // double r1 = vdiag.Length(); - double r1 = sqrt (3.0) * max3(vdiag.X(), vdiag.Y(), vdiag.Z()); - vdiag = Vec3d (r1, r1, r1); + double r1 = sqrt (3.0) * max3(vdiag(0), vdiag(1), vdiag(2)); + vdiag = Vec<3> (r1, r1, r1); //double r2; - Point3d pmin2 = pmin - 8 * vdiag; - Point3d pmax2 = pmax + 8 * vdiag; + Point<3> pmin2 = pmin - 8 * vdiag; + Point<3> pmax2 = pmax + 8 * vdiag; - Point3d cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2); - cp2.X() = pmin2.X(); - cp3.Y() = pmin2.Y(); - cp4.Z() = pmin2.Z(); + Point<3> cp1(pmin2), cp2(pmax2), cp3(pmax2), cp4(pmax2); + cp2(0) = pmin2(0); + cp3(1) = pmin2(1); + cp4(2) = pmin2(2); - - - int np = mesh.GetNP(); + size_t np = mesh.GetNP(); startel[0] = mesh.AddPoint (cp1); startel[1] = mesh.AddPoint (cp2); @@ -655,24 +564,21 @@ namespace netgen // flag points to use for Delaunay: BitArrayChar usep(np); usep.Clear(); - for (int i = 1; i <= adfront->GetNF(); i++) - { - const MiniElement2d & face = adfront->GetFace(i); - for (int j = 0; j < face.GetNP(); j++) - usep.Set (face[j]); - } - for (int i = oldnp + PointIndex::BASE; + for (auto & face : adfront->Faces()) + for (PointIndex pi : face.Face().PNums()) + usep.Set (pi); + + for (size_t i = oldnp + PointIndex::BASE; i < np + PointIndex::BASE; i++) usep.Set (i); - for (int i = 0; i < mesh.LockedPoints().Size(); i++) - usep.Set (mesh.LockedPoints()[i]); - + for (PointIndex pi : mesh.LockedPoints()) + usep.Set (pi); + Array freelist; - int cntp = 0; MeshNB meshnb (tempels, mesh.GetNP() + 5); @@ -689,13 +595,12 @@ namespace netgen list.AddElement (1); Array connected, treesearch; - - tpmin = tpmax = mesh.Point(startel[0]); - for (int k = 1; k < 4; k++) - { - tpmin.SetToMin (mesh.Point (startel[k])); - tpmax.SetToMax (mesh.Point (startel[k])); - } + Box<3> tbox(Box<3>::EMPTY_BOX); + for (size_t k = 0; k < 4; k++) + tbox.Add (mesh.Point(startel[k])); + Point<3> tpmin = tbox.PMin(); + Point<3> tpmax = tbox.PMax(); + tpmax = tpmax + 0.01 * (tpmax - tpmin); tettree.Insert (tpmin, tpmax, 1); @@ -727,6 +632,7 @@ namespace netgen for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++) mixed[pi] = PointIndex ( (prim * pi) % np + PointIndex::BASE ); + tloop.Start(); for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End()-4; pi++) { if (pi % 1000 == 0) @@ -755,7 +661,8 @@ namespace netgen connected, treesearch, freelist, list, insphere, closesphere); } - + tloop.Stop(); + for (int i = tempels.Size(); i >= 1; i--) if (tempels.Get(i)[0] <= 0) tempels.DeleteElement (i); @@ -783,6 +690,8 @@ namespace netgen void Meshing3 :: Delaunay (Mesh & mesh, int domainnr, const MeshingParameters & mp) { + static Timer t("Meshing3::Delaunay"); RegionTimer reg(t); + int np, ne; PrintMessage (1, "Delaunay meshing"); @@ -810,26 +719,27 @@ namespace netgen // improve delaunay - mesh by swapping !!!! Mesh tempmesh; - for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) - tempmesh.AddPoint (mesh[pi]); + + for (auto & meshpoint : mesh.Points()) + tempmesh.AddPoint (meshpoint); - for (int i = 1; i <= tempels.Size(); i++) + for (auto & tempel : tempels) { Element el(4); for (int j = 0; j < 4; j++) - el[j] = tempels.Elem(i)[j]; + el[j] = tempel[j]; el.SetIndex (1); - const Point3d & lp1 = mesh.Point (el[0]); - const Point3d & lp2 = mesh.Point (el[1]); - const Point3d & lp3 = mesh.Point (el[2]); - const Point3d & lp4 = mesh.Point (el[3]); - Vec3d v1(lp1, lp2); - Vec3d v2(lp1, lp3); - Vec3d v3(lp1, lp4); + const Point<3> & lp1 = mesh.Point (el[0]); + const Point<3> & lp2 = mesh.Point (el[1]); + const Point<3> & lp3 = mesh.Point (el[2]); + const Point<3> & lp4 = mesh.Point (el[3]); + Vec<3> v1 = lp2-lp1; + Vec<3> v2 = lp3-lp1; + Vec<3> v3 = lp4-lp1; - Vec3d n = Cross (v1, v2); + Vec<3> n = Cross (v1, v2); double vol = n * v3; if (vol > 0) swap (el[2], el[3]); @@ -858,7 +768,7 @@ namespace netgen { Element2d self(TRIG); self.SetIndex (1); - startel.GetFace1 (i, self); + startel.GetFace (i-1, self); tempmesh.AddSurfaceElement (self); } @@ -1232,6 +1142,7 @@ namespace netgen INDEX_3_HASHTABLE boundaryfaces(mesh.GetNOpenElements()/3+1); + /* for (int i = 1; i <= mesh.GetNOpenElements(); i++) { const Element2d & tri = mesh.OpenElement(i); @@ -1239,6 +1150,13 @@ namespace netgen i3.Sort(); boundaryfaces.PrepareSet (i3); } + */ + for (const Element2d & tri : mesh.OpenElements()) + { + INDEX_3 i3 (tri[0], tri[1], tri[2]); + i3.Sort(); + boundaryfaces.PrepareSet (i3); + } boundaryfaces.AllocateElements(); for (int i = 1; i <= mesh.GetNOpenElements(); i++) { @@ -1248,14 +1166,23 @@ namespace netgen boundaryfaces.Set (i3, 1); } + /* for (int i = 0; i < tempels.Size(); i++) for (int j = 0; j < 4; j++) tempels[i].NB(j) = 0; - + */ + for (auto & el : tempels) + for (int j = 0; j < 4; j++) + el.NB(j) = 0; + TABLE elsonpoint(mesh.GetNP()); + /* for (int i = 0; i < tempels.Size(); i++) { const DelaunayTet & el = tempels[i]; + */ + for (const DelaunayTet & el : tempels) + { INDEX_4 i4(el[0], el[1], el[2], el[3]); i4.Sort(); elsonpoint.IncSizePrepare (i4.I1()); @@ -1279,7 +1206,8 @@ namespace netgen INDEX_3_CLOSED_HASHTABLE faceht(100); Element2d hel(TRIG); - for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) + // for (PointIndex pi = mesh.Points().Begin(); pi < mesh.Points().End(); pi++) + for (PointIndex pi : mesh.Points().Range()) { faceht.SetSize (4 * elsonpoint[pi].Size()); for (int ii = 0; ii < elsonpoint[pi].Size(); ii++) @@ -1289,7 +1217,7 @@ namespace netgen for (int j = 1; j <= 4; j++) { - el.GetFace1 (j, hel); + el.GetFace (j-1, hel); hel.Invert(); hel.NormalizeNumbering(); @@ -1303,8 +1231,8 @@ namespace netgen { INDEX_2 i2 = faceht.Get(i3); - tempels.Elem(i).NB1(j) = i2.I1(); - tempels.Elem(i2.I1()).NB1(i2.I2()) = i; + tempels.Elem(i).NB(j-1) = i2.I1(); + tempels.Elem(i2.I1()).NB(i2.I2()-1) = i; } else { @@ -1464,7 +1392,7 @@ namespace netgen for (int j = 1; j <= 4; j++) { - INDEX_3 i3 = tempels.Get(ei).GetFace1(j); + INDEX_3 i3 = tempels.Get(ei).GetFace(j-1); /* Element2d face; tempels.Get(ei).GetFace(j, face); @@ -1474,8 +1402,8 @@ namespace netgen i3.Sort(); - if (tempels.Get(ei).NB1(j)) - elstack.Append (tempels.Get(ei).NB1(j)); + if (tempels.Get(ei).NB(j-1)) + elstack.Append (tempels.Get(ei).NB(j-1)); /* if (innerfaces.Used(i3)) diff --git a/libsrc/meshing/global.cpp b/libsrc/meshing/global.cpp index df6bcece..327727be 100644 --- a/libsrc/meshing/global.cpp +++ b/libsrc/meshing/global.cpp @@ -31,6 +31,9 @@ namespace netgen DLL_HEADER shared_ptr ng_geometry; // TraceGlobal glob2("global2"); + // global communicator for netgen + DLL_HEADER MPI_Comm ng_comm = MPI_COMM_WORLD; + weak_ptr global_mesh; void SetGlobalMesh (shared_ptr m) { @@ -50,10 +53,6 @@ namespace netgen string ngdir = "."; - // parallel netgen - int id = 0, ntasks = 1; - - void Ng_PrintDest(const char * s) { if (id == 0) diff --git a/libsrc/meshing/global.hpp b/libsrc/meshing/global.hpp index e6cb59de..c97815e3 100644 --- a/libsrc/meshing/global.hpp +++ b/libsrc/meshing/global.hpp @@ -59,6 +59,10 @@ namespace netgen DLL_HEADER extern weak_ptr global_mesh; DLL_HEADER void SetGlobalMesh (shared_ptr m); + + // global communicator for netgen (dummy if no MPI) + extern DLL_HEADER MPI_Comm ng_comm; + } #endif diff --git a/libsrc/meshing/hprefinement.cpp b/libsrc/meshing/hprefinement.cpp index a26945e6..d2752ba8 100644 --- a/libsrc/meshing/hprefinement.cpp +++ b/libsrc/meshing/hprefinement.cpp @@ -22,8 +22,9 @@ namespace netgen { pnums[i] = -1; param[i][0] = param[i][1] = param[i][2] = 0; - domin=-1; domout=-1; // he: } + domin=-1; domout=-1; // he: + levelx = 0; levely = 0; levelz = 0; } HPRefElement :: HPRefElement () @@ -31,46 +32,40 @@ namespace netgen Reset(); } - HPRefElement :: HPRefElement(Element & el) + HPRefElement :: HPRefElement(Element & el) : + np(el.GetNV()), index(el.GetIndex()), levelx(0), levely(0), levelz(0), type(HP_NONE), domin(-1), domout(-1) //domin,out for segements { //Reset(); - np = el.GetNV(); for (int i=0; igeom) { @@ -829,18 +806,27 @@ namespace netgen HPRefElement newel(el); newel.type = hprs->neweltypes[j]; - // newel.index = elements[i].index; // newel.coarse_elnr = elements[i].coarse_elnr; - newel.levelx = newel.levely = newel.levelz = newlevel; + if (newel.type == HP_SEGM || + newel.type == HP_TRIG || + newel.type == HP_QUAD || + newel.type == HP_TET || + newel.type == HP_PRISM || + newel.type == HP_HEX || + newel.type == HP_PYRAMID) + newel.levelx = newel.levely = newel.levelz = newlevel; + else + newel.levelx = newel.levely = newel.levelz = newlevel+1; + switch(hprsnew->geom) { - case HP_SEGM: newel.np=2; break; - case HP_QUAD: newel.np=4; break; - case HP_TRIG: newel.np=3; break; - case HP_HEX: newel.np=8; break; - case HP_PRISM: newel.np=6; break; - case HP_TET: newel.np=4; break; + case HP_SEGM: newel.np=2; break; + case HP_QUAD: newel.np=4; break; + case HP_TRIG: newel.np=3; break; + case HP_HEX: newel.np=8; break; + case HP_PRISM: newel.np=6; break; + case HP_TET: newel.np=4; break; case HP_PYRAMID: newel.np=5; break; default: throw NgException (string("hprefinement.cpp: illegal type")); @@ -868,7 +854,7 @@ namespace netgen if (j == 0) elements[i] = newel; // overwrite old element - else + else elements.Append (newel); j++; } @@ -1337,7 +1323,7 @@ namespace netgen Array & hpelements = *mesh.hpelements; InitHPElements(mesh,hpelements); - + Array nplevel; nplevel.Append (mesh.GetNP()); diff --git a/libsrc/meshing/improve3.cpp b/libsrc/meshing/improve3.cpp index 0dde9ba6..943c4197 100644 --- a/libsrc/meshing/improve3.cpp +++ b/libsrc/meshing/improve3.cpp @@ -22,6 +22,8 @@ namespace netgen void MeshOptimize3d :: CombineImprove (Mesh & mesh, OPTIMIZEGOAL goal) { + static Timer t("MeshOptimize3d::CombineImprove"); RegionTimer reg(t); + int np = mesh.GetNP(); int ne = mesh.GetNE(); @@ -274,6 +276,8 @@ void MeshOptimize3d :: CombineImprove (Mesh & mesh, void MeshOptimize3d :: SplitImprove (Mesh & mesh, OPTIMIZEGOAL goal) { + static Timer t("MeshOptimize3d::SplitImprove"); RegionTimer reg(t); + double bad1, bad2, badmax, badlimit; int cnt = 0; @@ -569,6 +573,8 @@ void MeshOptimize3d :: SplitImprove (Mesh & mesh, void MeshOptimize3d :: SwapImprove (Mesh & mesh, OPTIMIZEGOAL goal, const BitArray * working_elements) { + static Timer t("MeshOptimize3d::SwapImprove"); RegionTimer reg(t); + PointIndex pi3(0), pi4(0), pi5(0), pi6(0); int cnt = 0; @@ -2303,6 +2309,8 @@ void MeshOptimize3d :: SwapImproveSurface (Mesh & mesh, OPTIMIZEGOAL goal, void MeshOptimize3d :: SwapImprove2 (Mesh & mesh, OPTIMIZEGOAL goal) { + static Timer t("MeshOptimize3d::SwapImprove2"); RegionTimer reg(t); + PointIndex pi1(0), pi2(0), pi3(0), pi4(0), pi5(0); Element el21(TET), el22(TET), el31(TET), el32(TET), el33(TET); diff --git a/libsrc/meshing/meshclass.cpp b/libsrc/meshing/meshclass.cpp index 40a399eb..3ba4c913 100644 --- a/libsrc/meshing/meshclass.cpp +++ b/libsrc/meshing/meshclass.cpp @@ -4988,6 +4988,10 @@ namespace netgen bool build_searchtree, const bool allowindex) const { + // const double pointtol = 1e-12; + // netgen::Point<3> pmin = p - Vec<3> (pointtol, pointtol, pointtol); + // netgen::Point<3> pmax = p + Vec<3> (pointtol, pointtol, pointtol); + if (dimension == 2 || (dimension==3 && !GetNE() && GetNSE())) { int ne; diff --git a/libsrc/meshing/meshclass.hpp b/libsrc/meshing/meshclass.hpp index fe72de2a..d1f09b91 100644 --- a/libsrc/meshing/meshclass.hpp +++ b/libsrc/meshing/meshclass.hpp @@ -32,10 +32,8 @@ namespace netgen /// point coordinates T_POINTS points; -#ifdef PARALLEL - // The communicator for this mesh. (more or less dummy for now!) + // The communicator for this mesh. Just a dummy if compiled without MPI. MPI_Comm comm; -#endif /// line-segments at edges Array segments; @@ -132,8 +130,8 @@ namespace netgen /// mesh access semaphors. NgMutex majormutex; - SYMBOLTABLE< Array* > userdata_int; - SYMBOLTABLE< Array* > userdata_double; + SymbolTable< Array* > userdata_int; + SymbolTable< Array* > userdata_double; mutable Array< Point3d > pointcurves; @@ -223,7 +221,7 @@ namespace netgen DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer = 1); DLL_HEADER PointIndex AddPoint (const Point3d & p, int layer, POINTTYPE type); - int GetNP () const { return points.Size(); } + auto GetNP () const { return points.Size(); } // [[deprecated("Use Point(PointIndex) instead of int !")]] MeshPoint & Point(int i) { return points.Elem(i); } @@ -293,7 +291,7 @@ namespace netgen timestamp = NextTimeStamp(); } - int GetNSE () const { return surfelements.Size(); } + auto GetNSE () const { return surfelements.Size(); } // [[deprecated("Use SurfaceElement(SurfaceElementIndex) instead of int !")]] Element2d & SurfaceElement(int i) { return surfelements.Elem(i); } @@ -318,7 +316,7 @@ namespace netgen // write to pre-allocated container, thread-safe DLL_HEADER void SetVolumeElement (ElementIndex sei, const Element & el); - int GetNE () const { return volelements.Size(); } + auto GetNE () const { return volelements.Size(); } // [[deprecated("Use VolumeElement(ElementIndex) instead of int !")]] Element & VolumeElement(int i) { return volelements.Elem(i); } @@ -456,7 +454,8 @@ namespace netgen const Element2d & OpenElement(int i) const { return openelements.Get(i); } - + auto & OpenElements() const { return openelements; } + /// are also quads open elements bool HasOpenQuads () const; @@ -606,6 +605,9 @@ namespace netgen int AddEdgeDescriptor(const EdgeDescriptor & fd) { edgedecoding.Append(fd); return edgedecoding.Size() - 1; } + MPI_Comm GetCommunicator() const { return this->comm; } + void SetCommunicator(MPI_Comm acomm); + /// DLL_HEADER void SetMaterial (int domnr, const string & mat); /// @@ -858,7 +860,11 @@ namespace netgen void SendMesh ( ) const; // Mesh * mastermesh, Array & neloc) const; /// loads a mesh sent from master processor void ReceiveParallelMesh (); - +#else + void Distribute () {} + void SendRecvMesh () {} + void Distribute (Array & volume_weights, Array & surface_weights, + Array & segment_weights){ } #endif diff --git a/libsrc/meshing/meshfunc.cpp b/libsrc/meshing/meshfunc.cpp index d1629d62..e3fab9b8 100644 --- a/libsrc/meshing/meshfunc.cpp +++ b/libsrc/meshing/meshfunc.cpp @@ -14,6 +14,8 @@ namespace netgen // extern double teterrpow; MESHING3_RESULT MeshVolume (MeshingParameters & mp, Mesh& mesh3d) { + static Timer t("MeshVolume"); RegionTimer reg(t); + int oldne; int meshed; @@ -639,6 +641,8 @@ namespace netgen Mesh & mesh3d) // const CSGeometry * geometry) { + static Timer t("OptimizeVolume"); RegionTimer reg(t); + int i; PrintMessage (1, "Volume Optimization"); @@ -698,6 +702,8 @@ namespace netgen void RemoveIllegalElements (Mesh & mesh3d) { + static Timer t("RemoveIllegalElements"); RegionTimer reg(t); + int it = 10; int nillegal, oldn; diff --git a/libsrc/meshing/meshing3.cpp b/libsrc/meshing/meshing3.cpp index 16b49c53..e30454d1 100644 --- a/libsrc/meshing/meshing3.cpp +++ b/libsrc/meshing/meshing3.cpp @@ -168,12 +168,17 @@ int Meshing3 :: AddConnectedPair (const INDEX_2 & apair) MESHING3_RESULT Meshing3 :: GenerateMesh (Mesh & mesh, const MeshingParameters & mp) { - static int meshing3_timer = NgProfiler::CreateTimer ("Meshing3::GenerateMesh"); - static int meshing3_timer_a = NgProfiler::CreateTimer ("Meshing3::GenerateMesh a"); - static int meshing3_timer_b = NgProfiler::CreateTimer ("Meshing3::GenerateMesh b"); - static int meshing3_timer_c = NgProfiler::CreateTimer ("Meshing3::GenerateMesh c"); - static int meshing3_timer_d = NgProfiler::CreateTimer ("Meshing3::GenerateMesh d"); - NgProfiler::RegionTimer reg (meshing3_timer); + static Timer t("Meshing3::GenerateMesh"); RegionTimer reg(t); + static Timer meshing3_timer_a("Meshing3::GenerateMesh a", 2); + static Timer meshing3_timer_b("Meshing3::GenerateMesh b", 2); + static Timer meshing3_timer_c("Meshing3::GenerateMesh c", 1); + static Timer meshing3_timer_d("Meshing3::GenerateMesh d", 2); + // static int meshing3_timer = NgProfiler::CreateTimer ("Meshing3::GenerateMesh"); + // static int meshing3_timer_a = NgProfiler::CreateTimer ("Meshing3::GenerateMesh a"); + // static int meshing3_timer_b = NgProfiler::CreateTimer ("Meshing3::GenerateMesh b"); + // static int meshing3_timer_c = NgProfiler::CreateTimer ("Meshing3::GenerateMesh c"); + // static int meshing3_timer_d = NgProfiler::CreateTimer ("Meshing3::GenerateMesh d"); + // NgProfiler::RegionTimer reg (meshing3_timer); Array locpoints; // local points @@ -269,20 +274,16 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) } const MiniElement2d & bel = adfront->GetFace (baseelem); - const Point3d & p1 = adfront->GetPoint (bel[0]); - const Point3d & p2 = adfront->GetPoint (bel[1]); - const Point3d & p3 = adfront->GetPoint (bel[2]); + const Point<3> p1 = adfront->GetPoint (bel[0]); + const Point<3> p2 = adfront->GetPoint (bel[1]); + const Point<3> p3 = adfront->GetPoint (bel[2]); + - // (*testout) << endl << "base = " << bel << endl; - - - Point3d pmid = Center (p1, p2, p3); + Point<3> pmid = Center (p1, p2, p3); double his = (Dist (p1, p2) + Dist(p1, p3) + Dist(p2, p3)) / 3; - double hshould; - - hshould = mesh.GetH (pmid); - + double hshould = mesh.GetH (pmid); + if (adfront->GetFace (baseelem).GetNP() == 4) hshould = max2 (his, hshould); @@ -292,13 +293,13 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) double hinner = hmax * (1 + stat.qualclass); double houter = hmax * (1 + 2 * stat.qualclass); - NgProfiler::StartTimer (meshing3_timer_a); + meshing3_timer_a.Start(); stat.qualclass = adfront -> GetLocals (baseelem, locpoints, locfaces, pindex, findex, connectedpairs, houter, hinner, locfacesplit); - NgProfiler::StopTimer (meshing3_timer_a); + meshing3_timer_a.Stop(); // (*testout) << "locfaces = " << endl << locfaces << endl; @@ -320,9 +321,6 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) } - - - if (testmode) { (*testout) << "baseelem = " << baseelem << " qualclass = " << stat.qualclass << endl; @@ -479,7 +477,8 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) (*testout) << endl; } - NgProfiler::StartTimer (meshing3_timer_c); + // NgProfiler::StartTimer (meshing3_timer_c); + meshing3_timer_c.Start(); found = ApplyRules (plainpoints, allowpoint, locfaces, locfacesplit, connectedpairs, @@ -489,8 +488,8 @@ GenerateMesh (Mesh & mesh, const MeshingParameters & mp) if (found >= 0) impossible = 0; if (found < 0) found = 0; - - NgProfiler::StopTimer (meshing3_timer_c); + meshing3_timer_c.Stop(); + // NgProfiler::StopTimer (meshing3_timer_c); if (!found) loktestmode = 0; diff --git a/libsrc/meshing/meshtype.cpp b/libsrc/meshing/meshtype.cpp index baed182e..b28cc92d 100644 --- a/libsrc/meshing/meshtype.cpp +++ b/libsrc/meshing/meshtype.cpp @@ -1028,6 +1028,9 @@ namespace netgen case 6: typ = PRISM; break; case 8: typ = HEX; break; case 10: typ = TET10; break; + case 13: typ = PYRAMID13; break; + case 15: typ = PRISM15; break; + case 20: typ = HEX20; break; default: cerr << "Element::Element: unknown element with " << np << " points" << endl; } orderx = ordery = orderz = 1; @@ -1108,6 +1111,9 @@ namespace netgen case 6: typ = PRISM; break; case 8: typ = HEX; break; case 10: typ = TET10; break; + case 13: typ = PYRAMID13; break; + case 15: typ = PRISM15; break; + case 20: typ = HEX20; break; // default: break; cerr << "Element::SetNP unknown element with " << np << " points" << endl; @@ -1126,7 +1132,10 @@ namespace netgen case PRISM: np = 6; break; case HEX: np = 8; break; case TET10: np = 10; break; + case PYRAMID13: np = 13; break; case PRISM12: np = 12; break; + case PRISM15: np = 15; break; + case HEX20: np = 20; break; default: break; cerr << "Element::SetType unknown type " << int(typ) << endl; @@ -1955,7 +1964,27 @@ namespace netgen shape(4) = p(2); break; } - + case PYRAMID13: + { + T x = p(0); + T y = p(1); + T z = p(2); + z *= 1-1e-12; + shape[0] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (-2*x - z + 2)*(-2*y - z + 2))*(-0.5*x - 0.5*y - 0.5*z + 0.25); + shape[1] = (0.5*x - 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(-2*y - z + 2)); + shape[2] = (-z + z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*x + z)*(2*y + z))*(0.5*x + 0.5*y + 0.5*z - 0.75); + shape[3] = (-0.5*x + 0.5*y - 0.25)*(-z - z*(2*x + z - 1)*(2*y + z - 1)/(-z + 1) + (2*y + z)*(-2*x - z + 2)); + shape[4] = z*(2*z - 1); + shape[5] = 2*x*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shape[6] = 4*x*y*(-2*x - 2*z + 2)/(-2*z + 2); + shape[7] = 2*y*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-2*z + 2); + shape[8] = 4*x*y*(-2*y - 2*z + 2)/(-2*z + 2); + shape[9] = z*(-2*x - 2*z + 2)*(-2*y - 2*z + 2)/(-z + 1); + shape[10] = 2*x*z*(-2*y - 2*z + 2)/(-z + 1); + shape[11] = 4*x*y*z/(-z + 1); + shape[12] = 2*y*z*(-2*x - 2*z + 2)/(-z + 1); + break; + } case PRISM: { shape(0) = p(0) * (1-p(2)); @@ -1966,6 +1995,30 @@ namespace netgen shape(5) = (1-p(0)-p(1)) * p(2); break; } + case PRISM15: + { + T x = p(0); + T y = p(1); + T z = p(2); + T lam = 1-x-y; + T lamz = 1-z; + shape[0] = (2*x*x-x) * (2*lamz*lamz-lamz); + shape[1] = (2*y*y-y) * (2*lamz*lamz-lamz); + shape[2] = (2*lam*lam-lam) * (2*lamz*lamz-lamz); + shape[3] = (2*x*x-x) * (2*z*z-z); + shape[4] = (2*y*y-y) * (2*z*z-z); + shape[5] = (2*lam*lam-lam) * (2*z*z-z); + shape[6] = 4 * x * y * (2*lamz*lamz-lamz); + shape[7] = 4 * x * lam * (2*lamz*lamz-lamz); + shape[8] = 4 * y * lam * (2*lamz*lamz-lamz); + shape[9] = x * 4 * z * (1-z); + shape[10] = y * 4 * z * (1-z); + shape[11] = lam * 4 * z * (1-z); + shape[12] = 4 * x * y * (2*z*z-z); + shape[13] = 4 * x * lam * (2*z*z-z); + shape[14] = 4 * y * lam * (2*z*z-z); + break; + } case HEX: { shape(0) = (1-p(0))*(1-p(1))*(1-p(2)); @@ -1978,6 +2031,43 @@ namespace netgen shape(7) = (1-p(0))*( p(1))*( p(2)); break; } + case HEX20: + { + T x = p(0); + T y = p(1); + T z = p(2); + shape[0] = (1-x)*(1-y)*(1-z); + shape[1] = x *(1-y)*(1-z); + shape[2] = x * y *(1-z); + shape[3] = (1-x)* y *(1-z); + shape[4] = (1-x)*(1-y)*(z); + shape[5] = x *(1-y)*(z); + shape[6] = x * y *(z); + shape[7] = (1-x)* y *(z); + + T sigma[8]={(1-x)+(1-y)+(1-z),x+(1-y)+(1-z),x+y+(1-z),(1-x)+y+(1-z), + (1-x)+(1-y)+z,x+(1-y)+z,x+y+z,(1-x)+y+z}; + + static const int e[12][2] = + { + { 0, 1 }, { 2, 3 }, { 3, 0 }, { 1, 2 }, + { 4, 5 }, { 6, 7 }, { 7, 4 }, { 5, 6 }, + { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }, + }; + + for (int i = 0; i < 12; i++) + { + T lame = shape[e[i][0]]+shape[e[i][1]]; + T xi = sigma[e[i][1]]-sigma[e[i][0]]; + shape[8+i] = (1-xi*xi)*lame; + } + for (int i = 0; i < 12; i++) + { + shape[e[i][0]] -= 0.5 * shape[8+i]; + shape[e[i][1]] -= 0.5 * shape[8+i]; + } + break; + } default: throw NgException("Element :: GetNewShape not implemented for that element"); } diff --git a/libsrc/meshing/meshtype.hpp b/libsrc/meshing/meshtype.hpp index 5ded7330..df997e36 100644 --- a/libsrc/meshing/meshtype.hpp +++ b/libsrc/meshing/meshtype.hpp @@ -21,8 +21,8 @@ namespace netgen SEGMENT = 1, SEGMENT3 = 2, TRIG = 10, QUAD=11, TRIG6 = 12, QUAD6 = 13, QUAD8 = 14, TET = 20, TET10 = 21, - PYRAMID = 22, PRISM = 23, PRISM12 = 24, - HEX = 25 + PYRAMID = 22, PRISM = 23, PRISM12 = 24, PRISM15 = 27, PYRAMID13 = 28, + HEX = 25, HEX20 = 26 }; /* @@ -45,7 +45,7 @@ namespace netgen }; -#define ELEMENT_MAXPOINTS 12 +#define ELEMENT_MAXPOINTS 20 #define ELEMENT2D_MAXPOINTS 8 @@ -647,7 +647,7 @@ namespace netgen /// ELEMENT_TYPE typ; /// number of points (4..tet, 5..pyramid, 6..prism, 8..hex, 10..quad tet, 12..quad prism) - int np:5; + int np:6; /// class flagstruct { public: @@ -708,18 +708,21 @@ namespace netgen /// uint8_t GetNV() const { - __assume(typ >= TET && typ <= HEX); + __assume(typ >= TET && typ <= PYRAMID13); switch (typ) { case TET: case TET10: return 4; - case PRISM12: - case PRISM: + case PRISM12: + case PRISM15: + case PRISM: return 6; case PYRAMID: + case PYRAMID13: return 5; case HEX: + case HEX20: return 8; default: // not a 3D element #ifdef DEBUG @@ -794,10 +797,12 @@ namespace netgen { case TET: case TET10: return 4; - case PYRAMID: return 5; - case PRISM: + case PYRAMID: case PYRAMID13: return 5; + case PRISM: + case PRISM15: case PRISM12: return 5; - case HEX: return 6; + case HEX: case HEX20: + return 6; default: #ifdef DEBUG PrintSysError ("element3d::GetNFaces not implemented for typ", typ) diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp index f574f1f6..e9c040ef 100644 --- a/libsrc/meshing/parallelmesh.cpp +++ b/libsrc/meshing/parallelmesh.cpp @@ -35,9 +35,8 @@ namespace netgen void Mesh :: SendRecvMesh () { - int id, np; - MPI_Comm_rank(this->comm, &id); - MPI_Comm_size(this->comm, &np); + int id = MyMPI_GetId(GetCommunicator()); + int np = MyMPI_GetNTasks(GetCommunicator()); if (np == 1) { throw NgException("SendRecvMesh called, but only one rank in communicator!!"); @@ -73,16 +72,18 @@ namespace netgen { Array sendrequests; + MPI_Comm comm = GetCommunicator(); + int id = MyMPI_GetId(comm); + int ntasks = MyMPI_GetNTasks(comm); + int dim = GetDimension(); - MyMPI_Bcast(dim); + MyMPI_Bcast(dim, comm); const_cast(GetTopology()).Update(); PrintMessage ( 3, "Sending nr of elements"); - MPI_Comm comm = this->comm; - Array num_els_on_proc(ntasks); num_els_on_proc = 0; for (ElementIndex ei = 0; ei < GetNE(); ei++) @@ -285,7 +286,7 @@ namespace netgen for (int dest = 1; dest < ntasks; dest++) { FlatArray verts = verts_of_proc[dest]; - sendrequests.Append (MyMPI_ISend (verts, dest, MPI_TAG_MESH+1)); + sendrequests.Append (MyMPI_ISend (verts, dest, MPI_TAG_MESH+1, comm)); MPI_Datatype mptype = MeshPoint::MyGetMPIType(); @@ -301,7 +302,7 @@ namespace netgen MPI_Type_commit (&newtype); MPI_Request request; - MPI_Isend( &points[0], 1, newtype, dest, MPI_TAG_MESH+1, MPI_COMM_WORLD, &request); + MPI_Isend( &points[0], 1, newtype, dest, MPI_TAG_MESH+1, comm, &request); sendrequests.Append (request); } @@ -367,7 +368,7 @@ namespace netgen } Array req_per; for(int dest = 1; dest < ntasks; dest++) - req_per.Append(MyMPI_ISend(pp_data[dest], dest, MPI_TAG_MESH+1)); + req_per.Append(MyMPI_ISend(pp_data[dest], dest, MPI_TAG_MESH+1, comm)); MPI_Waitall(req_per.Size(), &req_per[0], MPI_STATUS_IGNORE); PrintMessage ( 3, "Sending Vertices - distprocs"); @@ -395,7 +396,7 @@ namespace netgen } for ( int dest = 1; dest < ntasks; dest ++ ) - sendrequests.Append (MyMPI_ISend (distpnums[dest], dest, MPI_TAG_MESH+1)); + sendrequests.Append (MyMPI_ISend (distpnums[dest], dest, MPI_TAG_MESH+1, comm)); @@ -425,7 +426,7 @@ namespace netgen } for (int dest = 1; dest < ntasks; dest ++ ) - sendrequests.Append (MyMPI_ISend (elementarrays[dest], dest, MPI_TAG_MESH+2)); + sendrequests.Append (MyMPI_ISend (elementarrays[dest], dest, MPI_TAG_MESH+2, comm)); PrintMessage ( 3, "Sending Face Descriptors" ); @@ -442,7 +443,7 @@ namespace netgen } for (int dest = 1; dest < ntasks; dest++) - sendrequests.Append (MyMPI_ISend (fddata, dest, MPI_TAG_MESH+3)); + sendrequests.Append (MyMPI_ISend (fddata, dest, MPI_TAG_MESH+3, comm)); /** Surface Elements **/ @@ -526,7 +527,7 @@ namespace netgen }); // distribute sel data for (int dest = 1; dest < ntasks; dest++) - sendrequests.Append (MyMPI_ISend(selbuf[dest], dest, MPI_TAG_MESH+4)); + sendrequests.Append (MyMPI_ISend(selbuf[dest], dest, MPI_TAG_MESH+4, comm)); /** Segments **/ @@ -577,8 +578,8 @@ namespace netgen osegs_both.Append(osegs1[l]); } for(int l = 0; l(1, &nbcs), k, MPI_TAG_MESH+6); - sendrequests[6*(k-1)+1] = MyMPI_ISend(bcname_sizes, k, MPI_TAG_MESH+6); - (void) MPI_Isend(compiled_bcnames, tot_bcsize, MPI_CHAR, k, MPI_TAG_MESH+6, MPI_COMM_WORLD, &sendrequests[6*(k-1)+2]); + sendrequests[6*(k-1)] = MyMPI_ISend(FlatArray(1, &nbcs), k, MPI_TAG_MESH+6, comm); + sendrequests[6*(k-1)+1] = MyMPI_ISend(bcname_sizes, k, MPI_TAG_MESH+6, comm); + (void) MPI_Isend(compiled_bcnames, tot_bcsize, MPI_CHAR, k, MPI_TAG_MESH+6, comm, &sendrequests[6*(k-1)+2]); } /** Send mat-names **/ @@ -719,9 +720,9 @@ namespace netgen for(int j=0;j(1, &nmats), k, MPI_TAG_MESH+6); - sendrequests[6*(k-1)+4] = MyMPI_ISend(mat_sizes, k, MPI_TAG_MESH+6); - (void) MPI_Isend(compiled_mats, tot_matsize, MPI_CHAR, k, MPI_TAG_MESH+6, MPI_COMM_WORLD, &sendrequests[6*(k-1)+5]); + sendrequests[6*(k-1)+3] = MyMPI_ISend(FlatArray(1, &nmats), k, MPI_TAG_MESH+6, comm); + sendrequests[6*(k-1)+4] = MyMPI_ISend(mat_sizes, k, MPI_TAG_MESH+6, comm); + (void) MPI_Isend(compiled_mats, tot_matsize, MPI_CHAR, k, MPI_TAG_MESH+6, comm, &sendrequests[6*(k-1)+5]); } /* now wait ... **/ @@ -731,7 +732,7 @@ namespace netgen PrintMessage( 3, "send mesh complete"); - MPI_Barrier(MPI_COMM_WORLD); + MPI_Barrier(comm); } @@ -750,14 +751,17 @@ namespace netgen int timer_sels = NgProfiler::CreateTimer ("Receive surface elements"); NgProfiler::RegionTimer reg(timer); + int id = MyMPI_GetId(GetCommunicator()); + int ntasks = MyMPI_GetNTasks(GetCommunicator()); + int dim; - MyMPI_Bcast(dim); + MyMPI_Bcast(dim, comm); SetDimension(dim); // Receive number of local elements int nelloc; MPI_Scatter (NULL, 0, MPI_INT, - &nelloc, 1, MPI_INT, 0, MPI_COMM_WORLD); + &nelloc, 1, MPI_INT, 0, comm); paralleltop -> SetNE (nelloc); // string st; @@ -766,8 +770,7 @@ namespace netgen NgProfiler::StartTimer (timer_pts); Array verts; - MyMPI_Recv (verts, 0, MPI_TAG_MESH+1); - + MyMPI_Recv (verts, 0, MPI_TAG_MESH+1, comm); int numvert = verts.Size(); paralleltop -> SetNV (numvert); @@ -787,11 +790,10 @@ namespace netgen MPI_Datatype mptype = MeshPoint::MyGetMPIType(); MPI_Status status; - MPI_Recv( &points[1], numvert, mptype, 0, MPI_TAG_MESH+1, MPI_COMM_WORLD, &status); + MPI_Recv( &points[1], numvert, mptype, 0, MPI_TAG_MESH+1, comm, &status); Array pp_data; - MyMPI_Recv(pp_data, 0, MPI_TAG_MESH+1); - + MyMPI_Recv(pp_data, 0, MPI_TAG_MESH+1, comm); int maxidentnr = pp_data[0]; auto & idents = GetIdentifications(); @@ -815,7 +817,7 @@ namespace netgen } Array dist_pnums; - MyMPI_Recv (dist_pnums, 0, MPI_TAG_MESH+1); + MyMPI_Recv (dist_pnums, 0, MPI_TAG_MESH+1, comm); for (int hi = 0; hi < dist_pnums.Size(); hi += 3) paralleltop -> @@ -828,7 +830,7 @@ namespace netgen Element el; Array elarray; - MyMPI_Recv (elarray, 0, MPI_TAG_MESH+2); + MyMPI_Recv (elarray, 0, MPI_TAG_MESH+2, comm); NgProfiler::RegionTimer reg(timer_els); @@ -848,7 +850,7 @@ namespace netgen { Array fddata; - MyMPI_Recv (fddata, 0, MPI_TAG_MESH+3); + MyMPI_Recv (fddata, 0, MPI_TAG_MESH+3, comm); for (int i = 0; i < fddata.Size(); i += 6) { int faceind = AddFaceDescriptor @@ -863,7 +865,7 @@ namespace netgen NgProfiler::RegionTimer reg(timer_sels); Array selbuf; - MyMPI_Recv ( selbuf, 0, MPI_TAG_MESH+4); + MyMPI_Recv ( selbuf, 0, MPI_TAG_MESH+4, comm); int ii = 0; int sel = 0; @@ -894,7 +896,7 @@ namespace netgen { Array segmbuf; - MyMPI_Recv ( segmbuf, 0, MPI_TAG_MESH+5); + MyMPI_Recv ( segmbuf, 0, MPI_TAG_MESH+5, comm); Segment seg; int globsegi; @@ -939,14 +941,14 @@ namespace netgen /** Recv bc-names **/ int nbcs; - MyMPI_Recv(nbcs, 0, MPI_TAG_MESH+6); + MyMPI_Recv(nbcs, 0, MPI_TAG_MESH+6, comm); Array bcs(nbcs); - MyMPI_Recv(bcs, 0, MPI_TAG_MESH+6); + MyMPI_Recv(bcs, 0, MPI_TAG_MESH+6, comm); int size_bc = 0; for(int k=0;k matsz(nmats); - MyMPI_Recv(matsz, 0, MPI_TAG_MESH+6); + MyMPI_Recv(matsz, 0, MPI_TAG_MESH+6, comm); int size_mats = 0; for(int k=0;k epart(ne), npart(nn); - idxtype nparts = ntasks-1; - idxtype edgecut; + idxtype nparts = MyMPI_GetNTasks(GetCommunicator())-1; - idxtype ncommon = 3; - METIS_PartMeshDual (&ne, &nn, &eptr[0], &eind[0], NULL, NULL, &ncommon, &nparts, - NULL, NULL, - &edgecut, &epart[0], &npart[0]); - - /* - METIS_PartMeshNodal (&ne, &nn, &eptr[0], &eind[0], NULL, NULL, &nparts, - NULL, NULL, - &edgecut, &epart[0], &npart[0]); - */ - PrintMessage (3, "metis complete"); - // cout << "done" << endl; - - for (int i = 0; i < GetNE(); i++) - VolumeElement(i+1).SetPartition(epart[i] + 1); - for (int i = 0; i < GetNSE(); i++) - SurfaceElement(i+1).SetPartition(epart[i+GetNE()] + 1); - for (int i = 0; i < GetNSeg(); i++) - LineSegment(i+1).SetPartition(epart[i+GetNE()+GetNSE()] + 1); + if (nparts == 1) + { + for (int i = 0; i < GetNE(); i++) + VolumeElement(i+1).SetPartition(1); + for (int i = 0; i < GetNSE(); i++) + SurfaceElement(i+1).SetPartition(1); + for (int i = 0; i < GetNSeg(); i++) + LineSegment(i+1).SetPartition(1); + } + else + + { + idxtype edgecut; + + idxtype ncommon = 3; + METIS_PartMeshDual (&ne, &nn, &eptr[0], &eind[0], NULL, NULL, &ncommon, &nparts, + NULL, NULL, + &edgecut, &epart[0], &npart[0]); + + /* + METIS_PartMeshNodal (&ne, &nn, &eptr[0], &eind[0], NULL, NULL, &nparts, + NULL, NULL, + &edgecut, &epart[0], &npart[0]); + */ + PrintMessage (3, "metis complete"); + // cout << "done" << endl; + + for (int i = 0; i < GetNE(); i++) + VolumeElement(i+1).SetPartition(epart[i] + 1); + for (int i = 0; i < GetNSE(); i++) + SurfaceElement(i+1).SetPartition(epart[i+GetNE()] + 1); + for (int i = 0; i < GetNSeg(); i++) + LineSegment(i+1).SetPartition(epart[i+GetNE()+GetNSE()] + 1); + } + + // surface elements attached to volume elements Array boundarypoints (GetNP()); boundarypoints = false; @@ -1120,12 +1140,9 @@ namespace netgen for (SurfaceElementIndex sei = 0; sei < GetNSE(); sei++) { const Element2d & el = (*this)[sei]; - cout << "surf-el " << sei << " verts: " << endl; for (int j = 0; j < el.GetNP(); j++) { - cout << el[j] << " "; f(el[j], sei); } - cout << endl; } }; auto loop_els_3d = [&](auto f) { @@ -1150,7 +1167,6 @@ namespace netgen if(boundarypoints[vertex]) cnt[vertex]++; }); - cout << "count: " << endl << cnt << endl; TABLE pnt2el(cnt); loop_els([&](auto vertex, int index) { @@ -1277,8 +1293,9 @@ namespace netgen // call it only for the master ! void Mesh :: Distribute (Array & volume_weights , Array & surface_weights, Array & segment_weights) { - MPI_Comm_size(MPI_COMM_WORLD, &ntasks); - MPI_Comm_rank(MPI_COMM_WORLD, &id); + MPI_Comm comm = GetCommunicator(); + int id = MyMPI_GetId(comm); + int ntasks = MyMPI_GetNTasks(comm); if (id != 0 || ntasks == 1 ) return; @@ -1368,7 +1385,20 @@ namespace netgen eptr.Append (eind.Size()); Array epart(ne), npart(nn); - idxtype nparts = ntasks-1; + idxtype nparts = MyMPI_GetNTasks(GetCommunicator())-1; + + if (nparts == 1) + { + for (int i = 0; i < GetNE(); i++) + VolumeElement(i+1).SetPartition(1); + for (int i = 0; i < GetNSE(); i++) + SurfaceElement(i+1).SetPartition(1); + for (int i = 0; i < GetNSeg(); i++) + LineSegment(i+1).SetPartition(1); + return; + } + + idxtype edgecut; diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp index 80c70121..346b5e25 100644 --- a/libsrc/meshing/paralleltop.cpp +++ b/libsrc/meshing/paralleltop.cpp @@ -24,6 +24,9 @@ namespace netgen void ParallelMeshTopology :: Reset () { *testout << "ParallelMeshTopology::Reset" << endl; + + int id = MyMPI_GetId(mesh.GetCommunicator()); + int ntasks = MyMPI_GetNTasks(mesh.GetCommunicator()); if ( ntasks == 1 ) return; @@ -206,6 +209,10 @@ namespace netgen // cout << "UpdateCoarseGrid" << endl; // if (is_updated) return; + MPI_Comm comm = mesh.GetCommunicator(); + int id = MyMPI_GetId(comm); + int ntasks = MyMPI_GetNTasks(comm); + Reset(); static int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid"); NgProfiler::RegionTimer reg(timer); @@ -222,14 +229,14 @@ namespace netgen // MPI_Barrier (MPI_COMM_WORLD); - MPI_Group MPI_GROUP_WORLD; + MPI_Group MPI_GROUP_comm; MPI_Group MPI_LocalGroup; MPI_Comm MPI_LocalComm; int process_ranks[] = { 0 }; - MPI_Comm_group (MPI_COMM_WORLD, &MPI_GROUP_WORLD); - MPI_Group_excl (MPI_GROUP_WORLD, 1, process_ranks, &MPI_LocalGroup); - MPI_Comm_create (MPI_COMM_WORLD, MPI_LocalGroup, &MPI_LocalComm); + MPI_Comm_group (comm, &MPI_GROUP_comm); + MPI_Group_excl (MPI_GROUP_comm, 1, process_ranks, &MPI_LocalGroup); + MPI_Comm_create (comm, MPI_LocalGroup, &MPI_LocalComm); if (id == 0) return; diff --git a/libsrc/meshing/python_mesh.cpp b/libsrc/meshing/python_mesh.cpp index 6360174a..34b00a0f 100644 --- a/libsrc/meshing/python_mesh.cpp +++ b/libsrc/meshing/python_mesh.cpp @@ -17,6 +17,22 @@ namespace netgen { extern bool netgen_executable_started; extern shared_ptr ng_geometry; +#ifdef PARALLEL + /** we need allreduce in python-wrapped communicators **/ + template + inline T MyMPI_AllReduceNG (T d, const MPI_Op & op = MPI_SUM, MPI_Comm comm = ng_comm) + { + T global_d; + MPI_Allreduce ( &d, &global_d, 1, MyGetMPIType(), op, comm); + return global_d; + } +#else + enum { MPI_SUM = 0, MPI_MIN = 1, MPI_MAX = 2 }; + typedef int MPI_Op; + template + inline T MyMPI_AllReduceNG (T d, const MPI_Op & op = MPI_SUM, MPI_Comm comm = ng_comm) + { return d; } +#endif } @@ -80,6 +96,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(py::self-py::self) .def(py::self+Vec<2>()) .def(py::self-Vec<2>()) + .def("__getitem__", [](Point<2>& self, int index) { return self[index]; }) ; py::class_> (m, "Point3d") @@ -88,6 +105,7 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(py::self-py::self) .def(py::self+Vec<3>()) .def(py::self-Vec<3>()) + .def("__getitem__", [](Point<2>& self, int index) { return self[index]; }) ; m.def ("Pnt", FunctionPointer @@ -111,6 +129,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(-py::self) .def(double()*py::self) .def("Norm", &Vec<2>::Length) + .def("__getitem__", [](Vec<2>& vec, int index) { return vec[index]; }) + .def("__len__", [](Vec<2>& /*unused*/) { return 2; }) ; py::class_> (m, "Vec3d") @@ -121,6 +141,8 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) .def(-py::self) .def(double()*py::self) .def("Norm", &Vec<3>::Length) + .def("__getitem__", [](Vec<3>& vec, int index) { return vec[index]; }) + .def("__len__", [](Vec<3>& /*unused*/) { return 3; }) ; m.def ("Vec", FunctionPointer @@ -234,37 +256,19 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) py::class_(m, "Element3D") .def(py::init([](int index, py::list vertices) { - Element * newel = nullptr; - if (py::len(vertices) == 4) - { - newel = new Element(TET); - for (int i = 0; i < 4; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 5) - { - newel = new Element(PYRAMID); - for (int i = 0; i < 5; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 6) - { - newel = new Element(PRISM); - for (int i = 0; i < 6; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else if (py::len(vertices) == 8) - { - newel = new Element(HEX); - for (int i = 0; i < 8; i++) - (*newel)[i] = py::extract(vertices[i])(); - newel->SetIndex(index); - } - else - throw NgException ("cannot create element"); + std::map types = {{4, TET}, + {5, PYRAMID}, + {6, PRISM}, + {8, HEX}, + {10, TET10}, + {13, PYRAMID13}, + {15, PRISM15}, + {20, HEX20}}; + int np = py::len(vertices); + auto newel = new Element(types[np]); + for(int i=0; i(vertices[i]); + newel->SetIndex(index); return newel; }), py::arg("index")=1,py::arg("vertices"), @@ -316,6 +320,13 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) (*newel)[i] = py::extract(vertices[i])(); newel->SetIndex(index); } + else if (py::len(vertices) == 8) + { + newel = new Element2d(QUAD8); + for(int i = 0; i<8; i++) + (*newel)[i] = py::extract(vertices[i])(); + newel->SetIndex(index); + } else throw NgException("Inconsistent number of vertices in Element2D"); return newel; @@ -483,18 +494,21 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) py::class_>(m, "Mesh") // .def(py::init<>("create empty mesh")) - .def(py::init( [] (int dim) + .def(py::init( [] (int dim, shared_ptr pycomm) { auto mesh = make_shared(); + mesh->SetCommunicator(pycomm!=nullptr ? pycomm->comm : netgen::ng_comm); mesh -> SetDimension(dim); SetGlobalMesh(mesh); // for visualization mesh -> SetGeometry (nullptr); return mesh; } ), - py::arg("dim")=3 + py::arg("dim")=3, py::arg("comm")=nullptr ) .def(NGSPickle()) - + .def_property_readonly("comm", [](const Mesh & amesh) + { return make_shared(amesh.GetCommunicator()); }, + "MPI-communicator the Mesh lives in") /* .def("__init__", [](Mesh *instance, int dim) @@ -506,17 +520,33 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) ) */ - .def("__str__", &ToString) .def_property_readonly("_timestamp", &Mesh::GetTimeStamp) + .def("Distribute", [](shared_ptr self, shared_ptr pycomm) { + MPI_Comm comm = pycomm!=nullptr ? pycomm->comm : self->GetCommunicator(); + self->SetCommunicator(comm); + if(MyMPI_GetNTasks(comm)==1) return self; + // if(MyMPI_GetNTasks(comm)==2) throw NgException("Sorry, cannot handle communicators with NP=2!"); + // cout << " rank " << MyMPI_GetId(comm) << " of " << MyMPI_GetNTasks(comm) << " called Distribute " << endl; + if(MyMPI_GetId(comm)==0) self->Distribute(); + else self->SendRecvMesh(); + return self; + }, py::arg("comm")=nullptr) + .def("Receive", [](shared_ptr pycomm) { + auto mesh = make_shared(); + mesh->SetCommunicator(pycomm->comm); + mesh->SendRecvMesh(); + return mesh; + }) .def("Load", FunctionPointer ([](Mesh & self, const string & filename) { istream * infile; + MPI_Comm comm = self.GetCommunicator(); + id = MyMPI_GetId(comm); + ntasks = MyMPI_GetNTasks(comm); + #ifdef PARALLEL - MPI_Comm_rank(MPI_COMM_WORLD, &id); - MPI_Comm_size(MPI_COMM_WORLD, &ntasks); - char* buf = nullptr; int strs = 0; if(id==0) { @@ -544,10 +574,10 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) } /** Scatter the geometry-string **/ - MPI_Bcast(&strs, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&strs, 1, MPI_INT, 0, comm); if(id!=0) buf = new char[strs]; - MPI_Bcast(buf, strs, MPI_CHAR, 0, MPI_COMM_WORLD); + MPI_Bcast(buf, strs, MPI_CHAR, 0, comm); if(id==0) delete infile; infile = new istringstream(string((const char*)buf, (size_t)strs)); @@ -902,6 +932,61 @@ DLL_HEADER void ExportNetgenMeshing(py::module &m) printmessage_importance = importance; return old; })); + + py::class_> (m, "MPI_Comm") + .def_property_readonly ("rank", &PyMPI_Comm::Rank) + .def_property_readonly ("size", &PyMPI_Comm::Size) + // .def_property_readonly ("rank", [](PyMPI_Comm & c) { cout << "rank for " << c.comm << endl; return c.Rank(); }) + // .def_property_readonly ("size", [](PyMPI_Comm & c) { cout << "size for " << c.comm << endl; return c.Size(); }) +#ifdef PARALLEL + .def("Barrier", [](PyMPI_Comm & c) { MPI_Barrier(c.comm); }) + .def("WTime", [](PyMPI_Comm & c) { return MPI_Wtime(); }) +#else + .def("Barrier", [](PyMPI_Comm & c) { }) + .def("WTime", [](PyMPI_Comm & c) { return -1.0; }) +#endif + .def("Sum", [](PyMPI_Comm & c, double x) { return MyMPI_AllReduceNG(x, MPI_SUM, c.comm); }) + .def("Min", [](PyMPI_Comm & c, double x) { return MyMPI_AllReduceNG(x, MPI_MIN, c.comm); }) + .def("Max", [](PyMPI_Comm & c, double x) { return MyMPI_AllReduceNG(x, MPI_MAX, c.comm); }) + .def("Sum", [](PyMPI_Comm & c, int x) { return MyMPI_AllReduceNG(x, MPI_SUM, c.comm); }) + .def("Min", [](PyMPI_Comm & c, int x) { return MyMPI_AllReduceNG(x, MPI_MIN, c.comm); }) + .def("Max", [](PyMPI_Comm & c, int x) { return MyMPI_AllReduceNG(x, MPI_MAX, c.comm); }) + .def("Sum", [](PyMPI_Comm & c, size_t x) { return MyMPI_AllReduceNG(x, MPI_SUM, c.comm); }) + .def("Min", [](PyMPI_Comm & c, size_t x) { return MyMPI_AllReduceNG(x, MPI_MIN, c.comm); }) + .def("Max", [](PyMPI_Comm & c, size_t x) { return MyMPI_AllReduceNG(x, MPI_MAX, c.comm); }) + .def("SubComm", [](PyMPI_Comm & c, std::vector proc_list) { + Array procs(proc_list.size()); + for (int i = 0; i < procs.Size(); i++) + procs[i] = proc_list[i]; + if (!procs.Contains(c.Rank())) + throw Exception("rank "+ToString(c.Rank())+" not in subcomm"); + MPI_Comm subcomm = MyMPI_SubCommunicator(c.comm, procs); + return make_shared(subcomm, true); + + /* + Array procs; + if (py::extract (proc_list).check()) { + py::list pylist = py::extract (proc_list)(); + procs.SetSize(py::len(pyplist)); + for (int i = 0; i < py::len(pylist); i++) + procs[i] = py::extract(pylist[i])(); + } + else { + throw Exception("SubComm needs a list!"); + } + if(!procs.Size()) { + cout << "warning, tried to construct empty communicator, returning MPI_COMM_NULL" << endl; + return make_shared(MPI_COMM_NULL); + } + else if(procs.Size()==2) { + throw Exception("Sorry, NGSolve cannot handle NP=2."); + } + MPI_Comm subcomm = MyMPI_SubCommunicator(c.comm, procs); + return make_shared(subcomm, true); + */ + }, py::arg("procs")); + ; + } PYBIND11_MODULE(libmesh, m) { diff --git a/libsrc/meshing/secondorder.cpp b/libsrc/meshing/secondorder.cpp index 861f99ed..49f91662 100644 --- a/libsrc/meshing/secondorder.cpp +++ b/libsrc/meshing/secondorder.cpp @@ -223,6 +223,45 @@ namespace netgen { 3, 4, 10 }, { 4, 5, 11 }, }; + static int betw_prism15[9][3] = + { + { 0, 1, 6 }, + { 0, 2, 7 }, + { 1, 2, 8 }, + { 0, 3, 9 }, + { 1, 4, 10 }, + { 2, 5, 11 }, + { 3, 4, 12 }, + { 3, 5, 13 }, + { 4, 5, 14 } + }; + static int betw_pyramid[8][3] = + { + { 0, 1, 5 }, + { 3, 2, 6 }, + { 3, 0, 7 }, + { 1, 2, 8 }, + { 0, 4, 9 }, + { 1, 4, 10 }, + { 2, 4, 11 }, + { 3, 4, 12 } + }; + static int betw_hex[12][3] = + { + { 0, 1, 8 }, + { 2, 3, 9 }, + { 3, 0, 10 }, + { 1, 2, 11 }, + { 4, 5, 12 }, + { 6, 7, 13 }, + { 7, 4, 14 }, + { 5, 6, 15 }, + { 0, 4, 16 }, + { 1, 5, 17 }, + { 2, 6, 18 }, + { 3, 7, 19 }, + }; + int (*betw)[3] = NULL; switch (el.GetType()) @@ -243,6 +282,29 @@ namespace netgen onp = 6; break; } + case PRISM15: + { + betw = betw_prism15; + newel.SetType(PRISM15); + onp = 6; + break; + } + case PYRAMID: + case PYRAMID13: + { + betw = betw_pyramid; + newel.SetType(PYRAMID13); + onp = 5; + break; + } + case HEX: + case HEX20: + { + betw = betw_hex; + newel.SetType (HEX20); + onp = 8; + break; + } default: PrintSysError ("MakeSecondOrder, illegal vol type ", int(el.GetType())); } diff --git a/libsrc/meshing/smoothing3.cpp b/libsrc/meshing/smoothing3.cpp index e16799dd..1353e156 100644 --- a/libsrc/meshing/smoothing3.cpp +++ b/libsrc/meshing/smoothing3.cpp @@ -1353,6 +1353,8 @@ void Mesh :: ImproveMesh (const CSG eometry & geometry, OPTIMIZEGOAL goal) void Mesh :: ImproveMesh (const MeshingParameters & mp, OPTIMIZEGOAL goal) { + static Timer t("Mesh::ImproveMesh"); RegionTimer reg(t); + int typ = 1; (*testout) << "Improve Mesh" << "\n"; diff --git a/libsrc/meshing/topology.cpp b/libsrc/meshing/topology.cpp index cddc731f..dc57c347 100644 --- a/libsrc/meshing/topology.cpp +++ b/libsrc/meshing/topology.cpp @@ -543,7 +543,7 @@ namespace netgen cnt = 0; ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); @@ -581,7 +581,9 @@ namespace netgen // accumulate number of edges int ned = edge2vert.Size(); - for (auto v : mesh->Points().Range()) + + // for (size_t v = 0; v < mesh->GetNV(); v++) + for (size_t v : cnt.Range()) { auto hv = cnt[v]; cnt[v] = ned; @@ -595,7 +597,7 @@ namespace netgen // for (PointIndex v = PointIndex::BASE; v < nv+PointIndex::BASE; v++) ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_CLOSED_HASHTABLE v2eht(2*max_edge_on_vertex+10); @@ -719,7 +721,7 @@ namespace netgen // for (auto v : mesh.Points().Range()) NgProfiler::StartTimer (timer2b1); ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); @@ -754,7 +756,9 @@ namespace netgen // accumulate number of faces int nfa = oldnfa; - for (auto v : mesh->Points().Range()) + // for (auto v : Range(mesh->GetNV())) // Points().Range()) + // for (size_t v = 0; v < mesh->GetNV(); v++) + for (auto v : cnt.Range()) { auto hv = cnt[v]; cnt[v] = nfa; @@ -765,7 +769,7 @@ namespace netgen // for (auto v : mesh.Points().Range()) ParallelForRange - (tm, mesh->Points().Size(), + (tm, mesh->GetNV(), // Points().Size(), [&] (size_t begin, size_t end) { INDEX_3_CLOSED_HASHTABLE vert2face(2*max_face_on_vertex+10); diff --git a/libsrc/meshing/topology.hpp b/libsrc/meshing/topology.hpp index 831a1ca5..b4e646d2 100644 --- a/libsrc/meshing/topology.hpp +++ b/libsrc/meshing/topology.hpp @@ -212,13 +212,16 @@ inline short int MeshTopology :: GetNVertices (ELEMENT_TYPE et) return 4; case PYRAMID: + case PYRAMID13: return 5; case PRISM: case PRISM12: + case PRISM15: return 6; case HEX: + case HEX20: return 8; // default: @@ -244,9 +247,11 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et) case QUAD: case QUAD6: - case QUAD8: return 4; + case QUAD8: + return 8; + case TET: return 4; case TET10: @@ -254,14 +259,21 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et) case PYRAMID: return 5; + case PYRAMID13: + return 13; case PRISM: - case PRISM12: return 6; + case PRISM12: + return 12; + case PRISM15: + return 15; case HEX: return 8; + case HEX20: + return 20; // default: // cerr << "Ng_ME_GetNVertices, illegal element type " << et << endl; } @@ -272,7 +284,7 @@ inline short int MeshTopology :: GetNPoints (ELEMENT_TYPE et) inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et) { - __assume(et >= SEGMENT && et <= HEX); + __assume(et >= SEGMENT && et <= PYRAMID13); switch (et) { case SEGMENT: @@ -293,13 +305,16 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et) return 6; case PYRAMID: + case PYRAMID13: return 8; case PRISM: case PRISM12: + case PRISM15: return 9; case HEX: + case HEX20: return 12; default: return 0; @@ -312,7 +327,7 @@ inline short int MeshTopology :: GetNEdges (ELEMENT_TYPE et) inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et) { - __assume(et >= SEGMENT && et <= HEX); + __assume(et >= SEGMENT && et <= PYRAMID13); switch (et) { case SEGMENT: @@ -333,13 +348,16 @@ inline short int MeshTopology :: GetNFaces (ELEMENT_TYPE et) return 4; case PYRAMID: + case PYRAMID13: return 5; case PRISM: case PRISM12: + case PRISM15: return 5; case HEX: + case HEX20: return 6; default: @@ -436,13 +454,16 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges1 (ELEMENT_TYPE et) return tet_edges; case PYRAMID: + case PYRAMID13: return pyramid_edges; case PRISM: case PRISM12: + case PRISM15: return prism_edges; case HEX: + case HEX20: return hex_edges; // default: // cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; @@ -513,7 +534,7 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) { 2, 6 }, { 3, 7 }, }; - + switch (et) { case SEGMENT: @@ -534,13 +555,16 @@ const ELEMENT_EDGE * MeshTopology :: GetEdges0 (ELEMENT_TYPE et) return tet_edges; case PYRAMID: + case PYRAMID13: return pyramid_edges; case PRISM: case PRISM12: + case PRISM15: return prism_edges; case HEX: + case HEX20: return hex_edges; // default: // cerr << "Ng_ME_GetEdges, illegal element type " << et << endl; @@ -618,15 +642,18 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces1 (ELEMENT_TYPE et) case PRISM: case PRISM12: + case PRISM15: return prism_faces; case PYRAMID: + case PYRAMID13: return pyramid_faces; case SEGMENT: case SEGMENT3: case HEX: + case HEX20: return hex_faces; // default: @@ -700,15 +727,18 @@ inline const ELEMENT_FACE * MeshTopology :: GetFaces0 (ELEMENT_TYPE et) case PRISM: case PRISM12: + case PRISM15: return prism_faces; case PYRAMID: + case PYRAMID13: return pyramid_faces; case SEGMENT: case SEGMENT3: case HEX: + case HEX20: return hex_faces; // default: diff --git a/libsrc/occ/CMakeLists.txt b/libsrc/occ/CMakeLists.txt index 831a3875..1f093ad8 100644 --- a/libsrc/occ/CMakeLists.txt +++ b/libsrc/occ/CMakeLists.txt @@ -8,11 +8,13 @@ if(USE_GUI) add_library(occvis ${NG_LIB_TYPE} vsocc.cpp) endif(USE_GUI) +target_link_libraries(occ PUBLIC ngcore) + if(NOT WIN32) - target_link_libraries( occ ${OCC_LIBRARIES} ${PYTHON_LIBRARIES}) + target_link_libraries( occ PUBLIC ${OCC_LIBRARIES} ${PYTHON_LIBRARIES}) install( TARGETS occ ${NG_INSTALL_DIR}) if (USE_GUI) - target_link_libraries( occvis occ ) + target_link_libraries( occvis PUBLIC occ ) install( TARGETS occvis ${NG_INSTALL_DIR}) endif(USE_GUI) endif(NOT WIN32) diff --git a/libsrc/occ/occmeshsurf.cpp b/libsrc/occ/occmeshsurf.cpp index a1f60db6..56fbf07b 100644 --- a/libsrc/occ/occmeshsurf.cpp +++ b/libsrc/occ/occmeshsurf.cpp @@ -673,13 +673,14 @@ namespace netgen if (surfi > 0) { - double u = gi1.u+secpoint*(gi2.u-gi1.u); double v = gi1.v+secpoint*(gi2.v-gi1.v); - - if (!geometry.FastProject (surfi, hnewp, u, v)) + + auto savept = hnewp; + if (!geometry.FastProject (surfi, hnewp, u, v) || Dist(hnewp, savept) > Dist(p1,p2)) { - // cout << "Fast projection to surface fails! Using OCC projection" << endl; + // cout << "Fast projection to surface fails! Using OCC projection" << endl; + hnewp = savept; geometry.Project (surfi, hnewp); } diff --git a/libsrc/stlgeom/CMakeLists.txt b/libsrc/stlgeom/CMakeLists.txt index 3620c496..8925c828 100644 --- a/libsrc/stlgeom/CMakeLists.txt +++ b/libsrc/stlgeom/CMakeLists.txt @@ -4,11 +4,13 @@ add_library(stl ${NG_LIB_TYPE} ) if(NOT WIN32) - target_link_libraries( stl mesh ${PYTHON_LIBRARIES}) - target_link_libraries( stl ${PYTHON_LIBRARIES}) + target_link_libraries( stl mesh ${PYTHON_LIBRARIES}) install( TARGETS stl ${NG_INSTALL_DIR}) endif(NOT WIN32) +target_link_libraries( stl ngcore ) + + if(USE_GUI) add_library(stlvis ${NG_LIB_TYPE} vsstl.cpp diff --git a/libsrc/stlgeom/meshstlsurface.cpp b/libsrc/stlgeom/meshstlsurface.cpp index 4e267193..02753fb8 100644 --- a/libsrc/stlgeom/meshstlsurface.cpp +++ b/libsrc/stlgeom/meshstlsurface.cpp @@ -598,7 +598,6 @@ void STLSurfaceMeshing1 (STLGeometry & geom, for (int fnr = 1; fnr <= mesh.GetNFD(); fnr++) { - if (fnr == 100) NgProfiler::ClearTimers(); if (!opensegsperface[fnr]) continue; if (multithread.terminate) return; diff --git a/libsrc/visualization/CMakeLists.txt b/libsrc/visualization/CMakeLists.txt index 8f803d92..4288862c 100644 --- a/libsrc/visualization/CMakeLists.txt +++ b/libsrc/visualization/CMakeLists.txt @@ -9,10 +9,8 @@ endif(USE_GUI) add_library(visual ${NG_LIB_TYPE} ${LIB_VISUAL_SOURCES}) -if(NOT WIN32) - target_link_libraries( visual ${PYTHON_LIBRARIES} ${MPI_CXX_LIBRARIES} ${OPENGL_LIBRARIES} ) - install( TARGETS visual ${NG_INSTALL_DIR}) -endif(NOT WIN32) +target_link_libraries( visual ngcore ${PYTHON_LIBRARIES} ${MPI_CXX_LIBRARIES} ${OPENGL_LIBRARIES} ) +install( TARGETS visual ${NG_INSTALL_DIR}) install(FILES meshdoc.hpp mvdraw.hpp diff --git a/libsrc/visualization/soldata.hpp b/libsrc/visualization/soldata.hpp index 56df5d65..74143fc0 100644 --- a/libsrc/visualization/soldata.hpp +++ b/libsrc/visualization/soldata.hpp @@ -1,7 +1,7 @@ #ifndef FILE_SOLDATA #define FILE_SOLDATA - +#include // for tAVX namespace netgen { @@ -103,10 +103,10 @@ namespace netgen #ifdef __SSE__ virtual bool GetMultiSurfValue (size_t selnr, size_t facetnr, size_t npts, - const tAVXd * xref, - const tAVXd * x, - const tAVXd * dxdxref, - tAVXd * values) + const ngsimd::tAVXd * xref, + const ngsimd::tAVXd * x, + const ngsimd::tAVXd * dxdxref, + ngsimd::tAVXd * values) { cerr << "GetMultiSurfVaue not overloaded for SIMD" << endl; return false; diff --git a/libsrc/visualization/vsmesh.cpp b/libsrc/visualization/vsmesh.cpp index d758c626..a6c5ef4c 100644 --- a/libsrc/visualization/vsmesh.cpp +++ b/libsrc/visualization/vsmesh.cpp @@ -2682,7 +2682,7 @@ namespace netgen for (ElementIndex ei = 0; ei < mesh->GetNE(); ei++) { const Element & el = (*mesh)[ei]; - if (el.GetType() == PYRAMID && !el.IsDeleted()) + if ((el.GetType() == PYRAMID || el.GetType() == PYRAMID13) && !el.IsDeleted()) { int i = ei + 1; diff --git a/libsrc/visualization/vssolution.cpp b/libsrc/visualization/vssolution.cpp index 2669c1c6..3b0b61cc 100644 --- a/libsrc/visualization/vssolution.cpp +++ b/libsrc/visualization/vssolution.cpp @@ -1373,7 +1373,7 @@ namespace netgen - if ( el.GetType() == QUAD || el.GetType() == QUAD6 ) + if ( el.GetType() == QUAD || el.GetType() == QUAD6 || el.GetType() == QUAD8 ) { bool curved = curv.IsSurfaceElementCurved (sei); @@ -2857,6 +2857,7 @@ namespace netgen } case PRISM: case PRISM12: + case PRISM15: { lami[0] = (1-lam3) * (1-lam1-lam2); lami[1] = (1-lam3) * lam1; @@ -2907,6 +2908,7 @@ namespace netgen } case PRISM: case PRISM12: + case PRISM15: { lami[0] = (1-lam3) * (1-lam1-lam2); lami[1] = (1-lam3) * lam1; @@ -2918,6 +2920,7 @@ namespace netgen break; } case PYRAMID: + case PYRAMID13: { if (lam3 > 1-1e-5) { @@ -3482,6 +3485,7 @@ namespace netgen case QUAD: case QUAD6: + case QUAD8: lami[0] = (1-lam1)*(1-lam2); lami[1] = lam1 * (1-lam2); lami[2] = lam1 * lam2; @@ -3723,6 +3727,7 @@ namespace netgen case QUAD: case QUAD6: + case QUAD8: lami[0] = (1-lam1)*(1-lam2); lami[1] = lam1 * (1-lam2); lami[2] = lam1 * lam2; @@ -4017,7 +4022,7 @@ namespace netgen if(vispar.donotclipdomain > 0 && vispar.donotclipdomain == (*mesh)[ei].GetIndex()) continue; ELEMENT_TYPE type = (*mesh)[ei].GetType(); - if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID) + if (type == HEX || type == PRISM || type == TET || type == TET10 || type == PYRAMID || type == PYRAMID13 || type == PRISM15 || type == HEX20) { const Element & el = (*mesh)[ei]; @@ -4076,6 +4081,8 @@ namespace netgen switch (type) { case PRISM: + case PRISM12: + case PRISM15: if (ix+iy <= n) { ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n); @@ -4086,9 +4093,11 @@ namespace netgen compress[ii] = -1; break; case HEX: + case HEX20: ploc = Point<3> (double(ix) / n, double(iy) / n, double(iz) / n); break; case PYRAMID: + case PYRAMID13: ploc = Point<3> (double(ix) / n * (1-double(iz)/n), double(iy) / n * (1-double(iz)/n), double(iz)/n); @@ -4102,7 +4111,7 @@ namespace netgen locgrid[compress[ii]] = ploc; } - if (type != TET && type != TET10 && type != PRISM) cnt_valid = n3; + if (type != TET && type != TET10 && type != PRISM && type != PRISM12 && type != PRISM15) cnt_valid = n3; locgrid.SetSize(cnt_valid); diff --git a/ng/ngpkg.cpp b/ng/ngpkg.cpp index a80b0ea6..5b67637a 100644 --- a/ng/ngpkg.cpp +++ b/ng/ngpkg.cpp @@ -1840,9 +1840,9 @@ namespace netgen - SYMBOLTABLE & GetVisualizationScenes () + SymbolTable & GetVisualizationScenes () { - static SYMBOLTABLE vss; + static SymbolTable vss; return vss; } @@ -1860,7 +1860,7 @@ namespace netgen vs = &vscross; if (GetVisualizationScenes().Used(vismode)) { - vs = GetVisualizationScenes().Get(vismode); + vs = GetVisualizationScenes()[vismode]; } else if (vismode) { diff --git a/nglib/CMakeLists.txt b/nglib/CMakeLists.txt index ff8bd1d5..7253230a 100644 --- a/nglib/CMakeLists.txt +++ b/nglib/CMakeLists.txt @@ -8,9 +8,6 @@ if(WIN32) $ $ $ - $ - $ - $ $ $ @@ -32,6 +29,7 @@ if(NOT WIN32) endif(USE_GUI) endif(NOT WIN32) +# target_link_libraries(nglib PRIVATE gen la gprim PUBLIC ngcore) target_link_libraries(nglib PUBLIC ngcore) target_link_libraries( nglib PRIVATE ${OCC_LIBRARIES} ${MPI_CXX_LIBRARIES} ${OPENGL_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT} ${X11_Xmu_LIB} ${JPEG_LIBRARIES} ${MKL_LIBRARIES} ${ZLIB_LIBRARIES} ${OCC_LIBRARIES} ) diff --git a/nglib/nglib.cpp b/nglib/nglib.cpp index 5fea184c..54fff096 100644 --- a/nglib/nglib.cpp +++ b/nglib/nglib.cpp @@ -1211,14 +1211,14 @@ namespace netgen } */ - + /* #ifndef WIN32 void ResetTime () { ; } #endif - + */ void MyBeep (int i) diff --git a/python/geom2d.py b/python/geom2d.py index 88856964..34b37870 100644 --- a/python/geom2d.py +++ b/python/geom2d.py @@ -39,11 +39,107 @@ def MakeCircle (geo, c, r, **args): geo.Append( ["spline3", pts[p1], pts[p2], pts[p3]], **args) + +def CreatePML(geo, pml_size, tol=1e-12): + """Create a pml layer around the geometry. This function works only on convex geometries and +the highest existing domain number must be named by using the function geo.SetMaterial(domnr, name). +Points in the geometry are assumed to be the same if (pt1 - pt2).Norm() < tol. +Returned is a dict with information to create the pml layer: + normals: A dict from the names of the linear pml domains to the normal vectors pointing inside the pml.""" + + def Start(spline): + if spline.rightdom == 0: + return spline.StartPoint() + return spline.EndPoint() + def End(spline): + if spline.rightdom == 0: + return spline.EndPoint() + return spline.StartPoint() + + splines = [] + for i in range(geo.GetNSplines()): + splines.append(geo.GetSpline(i)) + border = [] + is_closed = False + current_endpoint = None + while not is_closed: + for spline in splines: + if spline.leftdom == 0 or spline.rightdom == 0: + if current_endpoint is not None: + if (Start(spline)-current_endpoint).Norm() < tol: + border.append(spline) + current_endpoint = End(spline) + if (current_endpoint - startpoint).Norm() < tol: + is_closed = True + break + else: + startpoint = Start(spline) + current_endpoint = End(spline) + border.append(spline) + break + else: + raise Exception("Couldn't find closed spline around domain") + endpointindex_map = [] + for spline in border: + pnt = End(spline) + for i in range(geo.GetNPoints()): + if (pnt - geo.GetPoint(i)).Norm() < tol: + endpointindex_map.append(i) + break + else: + raise Exception("Couldn't find endpoint of spline in geometry") + start_ndoms = ndoms = geo.GetNDomains() + 1 + new_spline_domains = [] + normals = {} + for i, spline in enumerate(border): + if i == 0: + global_start = Start(spline) + pml_size * spline.GetNormal(0) + global_start_pnt = current_start = geo.AppendPoint(global_start[0], global_start[1]) + next_spline = border[(i+1)%len(border)] + new_end = End(spline) + pml_size * spline.GetNormal(1) + spline_name = geo.GetBCName(spline.bc) + if (new_end - global_start).Norm() < tol: + new_spline_domains.append(ndoms) + geo.Append(["line", current_start, global_start_pnt], bc="outer_" + spline_name, leftdomain = ndoms) + geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms) + geo.SetMaterial(ndoms, "pml_" + spline_name) + normals["pml_" + spline_name] = spline.GetNormal(0) + ndoms += 1 + break + end = geo.AppendPoint(new_end[0], new_end[1]) + new_spline_domains.append(ndoms) + geo.Append(["line", current_start, end], bc="outer_" + spline_name, leftdomain = ndoms) + geo.Append(["line", end, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1) + geo.SetMaterial(ndoms, "pml_" + spline_name) + normals["pml_" + spline_name] = spline.GetNormal(0) + ndoms += 1 + new_start = Start(next_spline) + pml_size * next_spline.GetNormal(0) + if (new_start - global_start).Norm() < tol: + geo.Append(["line", end, global_start_pnt], bc="outer", leftdomain = ndoms) + geo.Append(["line", global_start_pnt, endpointindex_map[i]], leftdomain=ndoms, rightdomain=start_ndoms) + geo.SetMaterial(ndoms, "pml_corner") + ndoms += 1 + break + if (new_end - new_start).Norm() < tol: + current_start = end + else: + current_start = geo.AppendPoint(new_start[0], new_start[1]) + geo.Append(["line", end, current_start], bc="outer", leftdomain = ndoms) + geo.Append(["line", current_start, endpointindex_map[i]], leftdomain=ndoms, rightdomain=ndoms+1) + geo.SetMaterial(ndoms, "pml_corner") + ndoms += 1 + for spline, domnr in zip(border, new_spline_domains): + if spline.leftdom == 0: + spline.leftdom = domnr + else: + spline.rightdom = domnr + return {"normals" : normals} + 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.AddSegment = lambda *args, **kwargs : SplineGeometry.Append(*args, **kwargs) SplineGeometry.AddPoint = lambda *args, **kwargs : SplineGeometry.AppendPoint(*args, **kwargs) - +SplineGeometry.CreatePML = CreatePML __all__ = ['SplineGeometry', 'unit_square'] diff --git a/python/read_gmsh.py b/python/read_gmsh.py index 0617732d..23a7ec52 100644 --- a/python/read_gmsh.py +++ b/python/read_gmsh.py @@ -1,9 +1,10 @@ from netgen.meshing import * - def ReadGmsh(filename): + if not filename.endswith(".msh"): + filename += ".msh" meshdim = 1 - with open(filename + ".msh", 'r') as f: + with open(filename, 'r') as f: while f.readline().split()[0] != "$Elements": pass nelem = int(f.readline()) @@ -16,15 +17,58 @@ def ReadGmsh(filename): meshdim = 3 break - f = open(filename + ".msh", 'r') + f = open(filename, 'r') mesh = Mesh(dim=meshdim) pointmap = {} facedescriptormap = {} - namemap = {} + namemap = { 0 : "default" } materialmap = {} bbcmap = {} + segm = 1 + trig = 2 + quad = 3 + tet = 4 + hex = 5 + prism = 6 + pyramid = 7 + segm3 = 8 # 2nd order line + trig6 = 9 # 2nd order trig + tet10 = 11 # 2nd order tet + point = 15 + quad8 = 16 # 2nd order quad + hex20 = 17 # 2nd order hex + prism15 = 18 # 2nd order prism + pyramid13 = 19 # 2nd order pyramid + segms = [segm, segm3] + trigs = [trig, trig6] + quads = [quad, quad8] + tets = [tet, tet10] + hexes = [hex, hex20] + prisms = [prism, prism15] + pyramids = [pyramid, pyramid13] + elem0d = [point] + elem1d = segms + elem2d = trigs + quads + elem3d = tets + hexes + prisms + pyramids + + num_nodes_map = { segm : 2, + trig : 3, + quad : 4, + tet : 4, + hex : 8, + prism : 6, + pyramid : 5, + segm3 : 3, + trig6 : 6, + tet10 : 10, + point : 1, + quad8 : 8, + hex20 : 20, + prism15 : 18, + pyramid13 : 19 } + while True: line = f.readline() if line == "": @@ -57,29 +101,14 @@ def ReadGmsh(filename): # the first tag is the physical group nr, the second tag is the group nr of the dim tags = [int(line[3 + k]) for k in range(numtags)] - if elmtype == 1: # 2-node line - num_nodes = 2 - elif elmtype == 2: # 3-node trig - num_nodes = 3 - elif elmtype == 3: # 4-node quad - num_nodes = 4 - elif elmtype == 4: # 4-node tet - num_nodes = 4 - elif elmtype == 5: # 8-node hex - num_nodes = 8 - elif elmtype == 6: # 6-node prism - num_nodes = 6 - elif elmtype == 7: # 5-node pyramid - num_nodes = 5 - elif elmtype == 15: # 1-node point - num_nodes = 1 - else: + if elmtype not in num_nodes_map: raise Exception("element type", elmtype, "not implemented") + num_nodes = num_nodes_map[elmtype] nodenums = line[3 + numtags:3 + numtags + num_nodes] nodenums2 = [pointmap[int(nn)] for nn in nodenums] - if elmtype == 1: + if elmtype in elem1d: if meshdim == 3: if tags[1] in bbcmap: index = bbcmap[tags[1]] @@ -117,7 +146,7 @@ def ReadGmsh(filename): mesh.Add(Element1D(index=index, vertices=nodenums2)) - if elmtype in [2, 3]: # 2d elements + if elmtype in elem2d: # 2d elements if meshdim == 3: if tags[1] in facedescriptormap.keys(): index = facedescriptormap[tags[1]] @@ -142,9 +171,17 @@ def ReadGmsh(filename): mesh.SetMaterial(index, "surf" + str(tags[1])) materialmap[tags[1]] = index - mesh.Add(Element2D(index, nodenums2)) + if elmtype in trigs: + ordering = [i for i in range(3)] + if elmtype == trig6: + ordering += [4,5,3] + if elmtype in quads: + ordering = [i for i in range(4)] + if elmtype == quad8: + ordering += [4, 6, 7, 5] + mesh.Add(Element2D(index, [nodenums2[i] for i in ordering])) - if elmtype in [4, 5, 6, 7]: # volume elements + if elmtype in elem3d: # volume elements if tags[1] in materialmap: index = materialmap[tags[1]] else: @@ -157,9 +194,22 @@ def ReadGmsh(filename): nodenums2 = [pointmap[int(nn)] for nn in nodenums] - if elmtype == 4: - mesh.Add(Element3D(index, [nodenums2[0], nodenums2[1], nodenums2[3], nodenums2[2]])) - elif elmtype in [5, 6, 7]: - mesh.Add(Element3D(index, nodenums2)) + if elmtype in tets: + ordering = [0,1,2,3] + if elmtype == tet10: + ordering += [4,6,7,5,9,8] + elif elmtype in hexes: + ordering = [0,1,5,4,3,2,6,7] + if elmtype == hex20: + ordering += [8,16,10,12,13,19,15,14,9,11,18,17] + elif elmtype in prisms: + ordering = [0,2,1,3,5,4] + if elmtype == prism15: + ordering += [7,6,9,8,11,10,13,12,14] + elif elmtype in pyramids: + ordering = [3,2,1,0,4] + if elmtype == pyramid13: + ordering += [10,5,6,8,12,11,9,7] + mesh.Add(Element3D(index, [nodenums2[i] for i in ordering])) return mesh diff --git a/tests/build_nospdlog.sh b/tests/build_nospdlog.sh new file mode 100644 index 00000000..2fedbef2 --- /dev/null +++ b/tests/build_nospdlog.sh @@ -0,0 +1,8 @@ +cd +mkdir -p build/netgen +cd build/netgen +cmake ../../src/netgen -DUSE_CCACHE=ON -DUSE_SPDLOG=OFF -DCMAKE_CXX_COMPILER=clang++ \ + -DCMAKE_C_COMPILER=clang +make -j12 +make install + diff --git a/tests/catch/CMakeLists.txt b/tests/catch/CMakeLists.txt index fc0ebe41..f0a25744 100644 --- a/tests/catch/CMakeLists.txt +++ b/tests/catch/CMakeLists.txt @@ -3,8 +3,7 @@ if(ENABLE_UNIT_TESTS) add_custom_target(unit_tests) # Build catch_main test object -message("netgen include dir = ${NETGEN_INCLUDE_DIR_ABSOLUTE} --------------------------------------") -include_directories(${CATCH_INCLUDE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/../../libsrc/include}) +include_directories(${CATCH_INCLUDE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/../../libsrc/include ${SPDLOG_INCLUDE_DIR}) add_library(catch_main STATIC main.cpp) set_target_properties(catch_main PROPERTIES CXX_STANDARD 17) add_dependencies(unit_tests catch_main) @@ -15,23 +14,19 @@ add_test(NAME unit_tests_built COMMAND ${CMAKE_COMMAND} --build . --target unit_ macro(add_unit_test name sources) add_executable(test_${name} ${sources} ) - if (WIN32) - target_link_libraries(test_${name} ngcore catch_main) - else(WIN32) - target_link_libraries(test_${name} ngcore catch_main) - endif(WIN32) + target_link_libraries(test_${name} ngcore catch_main) add_dependencies(unit_tests test_${name}) add_test(NAME unit_${name} COMMAND test_${name}) set_tests_properties(unit_${name} PROPERTIES DEPENDS unit_tests_built) + + if(ENABLE_CPP_CORE_GUIDELINES_CHECK) + set_target_properties(test_${name} PROPERTIES CXX_CLANG_TIDY "${DO_CLANG_TIDY}") + endif(ENABLE_CPP_CORE_GUIDELINES_CHECK) endmacro() add_unit_test(archive archive.cpp) +add_unit_test(symboltable symboltable.cpp) add_unit_test(version version.cpp) -if(ENABLE_CPP_CORE_GUIDELINES_CHECK) - set_target_properties(test_archive PROPERTIES CXX_CLANG_TIDY "${DO_CLANG_TIDY}") - set_target_properties(test_version PROPERTIES CXX_CLANG_TIDY "${DO_CLANG_TIDY}") -endif(ENABLE_CPP_CORE_GUIDELINES_CHECK) - endif(ENABLE_UNIT_TESTS) diff --git a/tests/catch/archive.cpp b/tests/catch/archive.cpp index b3f27d49..96c2087c 100644 --- a/tests/catch/archive.cpp +++ b/tests/catch/archive.cpp @@ -236,6 +236,34 @@ void testMultipleInheritance(Archive& in, Archive& out) } } +void testMap(Archive& in, Archive& out) +{ + map map1; + map1["netgen"] = "v6.2.1901"; + out & map1; + out.FlushBuffer(); + map map2; + in & map2; + CHECK(map2.size() == 1); + CHECK(map2["netgen"] == "v6.2.1901"); +} + +enum MyEnum + { + CASE1, + CASE2 + }; + +void testEnum(Archive& in, Archive& out) + { + MyEnum en = CASE2; + out & en; + out.FlushBuffer(); + MyEnum enin; + in & enin; + CHECK(enin == CASE2); + } + void testArchive(Archive& in, Archive& out) { SECTION("Empty String") @@ -285,11 +313,18 @@ void testArchive(Archive& in, Archive& out) { testNullPtr(in, out); } + SECTION("map") + { + testMap(in, out); + } + SECTION("enum") + { + testEnum(in, out); + } } TEST_CASE("BinaryArchive") { - SetLibraryVersion("netgen","v6.2.1811"); auto stream = make_shared(); BinaryOutArchive out(stream); BinaryInArchive in(stream); @@ -298,7 +333,6 @@ TEST_CASE("BinaryArchive") TEST_CASE("TextArchive") { - SetLibraryVersion("netgen","v6.2.1811"); auto stream = make_shared(); TextOutArchive out(stream); TextInArchive in(stream); diff --git a/tests/catch/symboltable.cpp b/tests/catch/symboltable.cpp new file mode 100644 index 00000000..5e6ecd2c --- /dev/null +++ b/tests/catch/symboltable.cpp @@ -0,0 +1,64 @@ + +#include "catch.hpp" +#include <../core/ngcore.hpp> +using namespace ngcore; +using namespace std; + +TEST_CASE("Symboltable") +{ + SymbolTable table; + CHECK_THROWS_AS(table["test"], RangeException); + table.Set("first", 1); + CHECK(table["first"] == 1); + table["first"] = 2; + CHECK(table["first"] == 2); + auto index = table.Index("first"); + CHECK(index == 0); + CHECK(table[index] == 2); + table[index] = 3; + CHECK(table["first"] == 3); +#ifndef NDEBUG + int a; + CHECK_THROWS_AS(a = table[5], RangeException); + CHECK_THROWS_AS(table.GetName(5), RangeException); +#endif + CHECK(table.Used("first")); + CHECK(!table.Used("second")); + SymbolTable another; + another.Set("first", 1); + another.Set("second", 2); + table.Update(another); + CHECK(table["first"] == 1); + CHECK(table["second"] == 2); + std::stringstream s; + s << table; + CHECK(s.str() == "first : 1\nsecond : 2\n"); + auto ss = std::make_shared(); + BinaryOutArchive ao(ss); + ao & table; + ao.FlushBuffer(); + BinaryInArchive ai(ss); + SymbolTable read; + ai & read; + for(size_t i = 0; i is special because of vector is special... + SymbolTable btable; + btable.Set("true", true); + btable.Set("false", false); + CHECK(btable[0]); + CHECK(!btable[1]); + CHECK(btable["true"]); + CHECK(!btable["false"]); + ao & btable; + ao.FlushBuffer(); + SymbolTable bread; + ai & bread; + CHECK(bread["true"]); + CHECK(!bread["false"]); +}