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

merge master

Browse Source
This commit is contained in:
Michael Neunteufel 2019-02-11 09:16:22 +01:00
commit e8f6ca5c1b
97 changed files with 4054 additions and 1338 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
.gitlab-ci.yml
View File

@ -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

17
CMakeLists.txt
View File

@ -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

3
cmake/SuperBuild.cmake
View File

@ -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

18
cmake/external_projects/spdlog.cmake Normal file
View File

@ -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)

2
cmake/external_projects/tcltk.cmake
View File

@ -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"

4
libsrc/core/.clang-tidy
View File

@ -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: '*'
WarningsAsErrors: '*'

39
libsrc/core/CMakeLists.txt
View File

@ -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 $<$<CONFIG:DEBUG>: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)

10
libsrc/core/archive.cpp
View File

@ -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<std::map<std::string, detail::ClassArchiveInfo>> type_register; // NOLINT
const detail::ClassArchiveInfo& Archive :: GetArchiveRegister(const std::string& classname)

255
libsrc/core/archive.hpp
View File

@ -3,31 +3,32 @@
#include <complex> // for complex
#include <cstring> // for size_t, strlen
#include <fstream> // for operator<<, ifstream, ofstream, basic...
#include <fstream> // for ifstream, ofstream
#include <functional> // for function
#include <map> // for map, _Rb_tree_iterator
#include <memory> // for __shared_ptr_access, __shared_ptr_acc...
#include <stdexcept> // for runtime_error
#include <string> // for string, operator+
#include <map> // for map
#include <memory> // for shared_ptr
#include <string> // for string
#include <type_traits> // for declval, enable_if, false_type, is_co...
#include <typeinfo> // for type_info
#include <utility> // for move, swap, pair
#include <vector> // 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 <pybind11/pybind11.h>
#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<typename T, typename ...Rest>
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<typename T, typename= typename std::enable_if<std::is_constructible<T>::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<void*, int> shared_ptr2nr, ptr2nr;
std::map<void*, int> shared_ptr2nr{}, ptr2nr{};
// vectors for storing the unarchived (shared) pointers
std::vector<std::shared_ptr<void>> nr2shared_ptr;
std::vector<void*> nr2ptr;
std::vector<std::shared_ptr<void>> nr2shared_ptr{};
std::vector<void*> nr2ptr{};
protected:
bool shallow_to_python = false;
std::map<std::string, VersionInfo> version_map = GetLibraryVersions();
std::shared_ptr<Logger> 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<typename T>
Archive& Shallow(T& val)
{
static_assert(detail::is_any_pointer<T>, "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<T>(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<typename T>
auto operator & (T& val) -> typename std::enable_if<std::is_enum<T>::value, Archive&>::type
{
int enumval;
if(Output())
enumval = int(val);
*this & enumval;
if(Input())
val = T(enumval);
return *this;
}
// vector<bool> 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<bool>& v)
{
logger->debug("In special archive for std::vector<bool>");
size_t size;
if(Output())
size = v.size();
(*this) & size;
if(Input())
{
v.resize(size);
bool b;
for(size_t i=0; i<size; i++)
{
(*this) & b;
v[i] = b;
}
}
else
{
for(bool b : v)
(*this) & b;
}
return *this;
}
template<typename T1, typename T2>
Archive& operator& (std::map<T1, T2>& 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<void*>(ptr.get()) )
if(reg_ptr != static_cast<void*>(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<void> 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<T>(other);
{
logger->debug("Shared pointer didn't need pointer casts");
ptr = std::static_pointer_cast<T>(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<void*>(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<void*>(p));
if(reg_ptr != static_cast<void*>(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<T>::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<void*>(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<T>();
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<void*>(dynamic_cast<B1*>(p))); }
catch(std::exception&)
catch(const Exception&)
{ return Caster<T, Brest...>::tryUpcast(ti, p); }
}
@ -520,7 +623,7 @@ namespace ngcore
return dynamic_cast<T*>(static_cast<B1*>(GetArchiveRegister(Demangle(typeid(B1).name())).
downcaster(ti, p)));
}
catch(std::exception&)
catch(const Exception&)
{
return Caster<T, Brest...>::tryDowncast(ti, p);
}
@ -536,7 +639,7 @@ namespace ngcore
{
static_assert(detail::all_of_tmpl<std::is_base_of<Bases,T>::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<T>()
: Archive::Caster<T, Bases...>::tryUpcast(ti, detail::constructIfPossible<T>()); };
@ -805,7 +908,7 @@ namespace ngcore
}
};
#ifdef NG_PYTHON
#ifdef NETGEN_PYTHON
template<typename ARCHIVE>
class PyArchive : public ARCHIVE
@ -813,7 +916,12 @@ namespace ngcore
private:
pybind11::list lst;
size_t index = 0;
std::map<std::string, VersionInfo> 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<std::stringstream>()),
@ -821,8 +929,30 @@ namespace ngcore
{
ARCHIVE::shallow_to_python = true;
if(Input())
stream = std::make_shared<std::stringstream>
(pybind11::cast<pybind11::bytes>(lst[pybind11::len(lst)-1]));
{
stream = std::make_shared<std::stringstream>
(pybind11::cast<pybind11::bytes>(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<std::stringstream>
(pybind11::cast<pybind11::bytes>(lst[pybind11::len(lst)-2]));
*this & version_map;
stream = std::make_shared<std::stringstream>
(pybind11::cast<pybind11::bytes>(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<std::stringstream>(stream)->str()));
stream = std::make_shared<std::stringstream>();
*this & GetLibraryVersions();
FlushBuffer();
lst.append(pybind11::bytes(std::static_pointer_cast<std::stringstream>(stream)->str()));
stream = std::make_shared<std::stringstream>();
logger->debug("Writeout version needed = {}", version_needed);
*this & version_needed;
FlushBuffer();
lst.append(pybind11::bytes(std::static_pointer_cast<std::stringstream>(stream)->str()));
return lst;
@ -843,27 +982,31 @@ namespace ngcore
};
template<typename T, typename T_ARCHIVE_OUT=BinaryOutArchive, typename T_ARCHIVE_IN=BinaryInArchive>
auto NGSPickle(bool printoutput=false)
auto NGSPickle()
{
return pybind11::pickle([printoutput](T* self)
return pybind11::pickle([](T* self)
{
PyArchive<T_ARCHIVE_OUT> 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<T_ARCHIVE_IN> ar(state[0]);
ar & val;
return val;
});
}
#endif // NG_PYTHON
#endif // NETGEN_PYTHON
} // namespace ngcore
#endif // NETGEN_CORE_ARCHIVE_HPP

90
libsrc/core/exception.hpp Normal file
View File

@ -0,0 +1,90 @@
#ifndef NETGEN_CORE_EXCEPTION_HPP
#define NETGEN_CORE_EXCEPTION_HPP
#include <sstream> // for stringstream
#include <stdexcept> // for exception
#include <string> // 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<typename T>
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

137
libsrc/core/logging.cpp Normal file
View File

@ -0,0 +1,137 @@
#include "logging.hpp"
#ifdef NETGEN_USE_SPDLOG
#include <spdlog/spdlog.h>
#include <spdlog/sinks/ansicolor_sink.h>
#include <spdlog/sinks/basic_file_sink.h>
#else // NETGEN_USE_SPDLOG
#include <iostream>
#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<std::shared_ptr<spdlog::sinks::sink>>& GetDefaultSinks()
{
static std::vector<std::shared_ptr<spdlog::sinks::sink>> sinks =
{ std::make_shared<spdlog::sinks::ansicolor_stdout_sink_mt>() };
return sinks;
}
std::shared_ptr<spdlog::logger> CreateDefaultLogger(const std::string& name)
{
auto& default_sinks = GetDefaultSinks();
auto logger = std::make_shared<spdlog::logger>(name, default_sinks.begin(), default_sinks.end());
spdlog::details::registry::instance().register_and_init(logger);
return logger;
}
} // namespace detail
std::shared_ptr<Logger> GetLogger(const std::string& name)
{
auto logger = spdlog::get(name);
if(!logger)
logger = detail::CreateDefaultLogger(name);
return std::make_shared<Logger>(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<spdlog::sinks::basic_file_sink_mt>(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<spdlog::sinks::ansicolor_stdout_sink_mt>();
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<Logger> GetLogger(const std::string& /*unused*/)
{
return std::make_shared<Logger>(std::make_shared<spdlog::logger>());
}
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

124
libsrc/core/logging.hpp Normal file
View File

@ -0,0 +1,124 @@
#ifndef NETGEN_CORE_LOGGING_HPP
#define NETGEN_CORE_LOGGING_HPP
#undef NETGEN_USE_SPDLOG
#include <iostream>
#include <memory>
#include <string>
#include <vector>
#include "exception.hpp"
#include "ngcore_api.hpp"
#include "utils.hpp"
#ifdef NETGEN_USE_SPDLOG
#include <spdlog/fmt/fmt.h>
#include <spdlog/fmt/ostr.h> // 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<spdlog::logger> logger;
Logger(std::shared_ptr<spdlog::logger> l) : logger(std::move(l)) {}
void NGCORE_API log( level::level_enum level, std::string && s);
#ifdef NETGEN_USE_SPDLOG
template<typename ... Args>
void log( level::level_enum level, const char* str, Args ... args)
{
log(level, fmt::format(str, args...));
}
#else // NETGEN_USE_SPDLOG
template<typename T>
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<typename T>
std::string log_helper(std::string s, const T &t)
{
return replace(s,t);
}
template<typename T, typename ... Args>
std::string log_helper( std::string s, const T &t, Args ... args)
{
return log_helper(replace(s,t), args...);
}
template<typename ... Args>
void log( level::level_enum level, const char* str, Args ... args)
{
log(level, log_helper(std::string(str), args...));
}
#endif // NETGEN_USE_SPDLOG
template<typename ... Args>
void trace( const char* str, Args ... args) { log(level::level_enum::trace, str, args...); }
template<typename ... Args>
void debug( const char* str, Args ... args) { log(level::level_enum::debug, str, args...); }
template<typename ... Args>
void info( const char* str, Args ... args) { log(level::level_enum::info, str, args...); }
template<typename ... Args>
void warn( const char* str, Args ... args) { log(level::level_enum::warn, str, args...); }
template<typename ... Args>
void error( const char* str, Args ... args) { log(level::level_enum::err, str, args...); }
template<typename ... Args>
void critical( const char* str, Args ... args) { log(level::level_enum::critical, str, args...); }
};
NGCORE_API std::shared_ptr<Logger> 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

4
libsrc/core/ngcore.hpp
View File

@ -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

28
libsrc/core/ngcore_api.hpp
View File

@ -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

667
libsrc/core/paje_trace.cpp Normal file
View File

@ -0,0 +1,667 @@
#include <algorithm>
#include <atomic>
#include <iostream>
#include <map>
#include <set>
#include <thread>
#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<size_t>(std::numeric_limits<int>::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<d)
r=1, g=6*x,b=0;
else if (x<2*d)
r=1.0-6*(x-d),g=1,b=0;
else if (x<3*d)
r=0, g=1,b=6*(x-2*d);
else if (x<4*d)
r=0, g=1-6*(x-3*d),b=1;
else if (x<5*d)
r=6*(x-4*d), g=0,b=1;
else
r=1, g=0,b=1-5*(x-d);
};
int alias_counter;
FILE * ctrace_stream;
TTimePoint start_time;
std::shared_ptr<Logger> logger = GetLogger("PajeTrace");
double ConvertTime(TTimePoint t) {
// return time in milliseconds as double
// return std::chrono::duration<double>(t-start_time).count()*1000.0;
// return std::chrono::duration<double>(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<PajeEvent> 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<std::string> shash;
size_t h = shash(name);
h ^= h>>32U;
h = static_cast<uint32_t>(h);
hue = h*1.0/std::numeric_limits<uint32_t>::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<int> container_nodes;
container_nodes.reserve(num_nodes);
for(int i=0; i<num_nodes; i++)
container_nodes.emplace_back( paje.CreateContainer( container_type_node, container_task_manager, "Node " + ToString(i)) );
std::vector <int> thread_aliases;
thread_aliases.reserve(nthreads);
if(trace_threads)
for (int i=0; i<nthreads; i++)
{
auto name = "Timer level " + ToString(i);
thread_aliases.emplace_back( paje.CreateContainer( container_type_thread, container_nodes[i*num_nodes/nthreads], name ) );
}
std::map<const std::type_info *, int> job_map;
std::map<const std::type_info *, int> 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<int> timer_ids;
std::map<int,int> 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<int> timer_container_aliases;
timer_container_aliases.resize(maxdepth);
for(int i=0; i<maxdepth; i++)
{
auto name = "Timer level " + ToString(i);
timer_container_aliases[i] = paje.CreateContainer( container_type_timer, container_task_manager, name );
}
timerdepth = 0;
for(auto & event : timer_events)
{
if(event.is_start)
paje.PushState( event.time, state_type_timer, timer_container_aliases[timerdepth++], timer_aliases[event.timer_id] );
else
paje.PopState( event.time, state_type_timer, timer_container_aliases[--timerdepth] );
}
for(auto & vtasks : tasks)
{
for (Task & t : vtasks) {
int value_id = t.id;
switch(t.id_type)
{
case Task::ID_JOB:
value_id = job_task_map[jobs[t.id-1].type];
if(trace_thread_counter)
{
paje.AddVariable( t.start_time, variable_type_active_threads, container_jobs, 1.0 );
paje.SubVariable( t.stop_time, variable_type_active_threads, container_jobs, 1.0 );
}
if(trace_threads)
{
paje.PushState( t.start_time, state_type_task, thread_aliases[t.thread_id], value_id, t.additional_value, true );
paje.PopState( t.stop_time, state_type_task, thread_aliases[t.thread_id] );
}
break;
case Task::ID_TIMER:
value_id = timer_aliases[t.id];
paje.PushState( t.start_time, state_type_timer, thread_aliases[t.thread_id], value_id, t.additional_value, true );
paje.PopState( t.stop_time, state_type_timer, thread_aliases[t.thread_id] );
break;
default:
paje.PushState( t.start_time, state_type_task, thread_aliases[t.thread_id], value_id, t.additional_value, false );
paje.PopState( t.stop_time, state_type_task, thread_aliases[t.thread_id] );
break;
}
}
}
// Merge link event
int nlinks = 0;
for( auto & l : links)
nlinks += l.size();
std::vector<ThreadLink> links_merged;
links_merged.reserve(nlinks);
std::vector<unsigned int> pos(nthreads);
int nlinks_merged = 0;
while(nlinks_merged < nlinks)
{
int minpos = -1;
TTimePoint mintime = -1;
for (int t = 0; t<nthreads; t++)
{
if(pos[t] < links[t].size() && (minpos==-1 || links[t][pos[t]].time < mintime))
{
minpos = t;
mintime = links[t][pos[t]].time;
}
}
links_merged.push_back( links[minpos][pos[minpos]] );
pos[minpos]++;
nlinks_merged++;
}
std::vector<ThreadLink> 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<started_links.size())
{
while(i<started_links.size() && started_links[i].key == l.key)
{
ThreadLink & sl = started_links[i];
// Avoid links on same thread
if(sl.thread_id != l.thread_id)
{
paje.StartLink( sl.time, link_type, container_nodes[sl.thread_id*num_nodes/nthreads], l.key, thread_aliases[sl.thread_id], l.key);
paje.EndLink( l.time, link_type, container_nodes[l.thread_id*num_nodes/nthreads], l.key, thread_aliases[l.thread_id], l.key);
}
started_links.erase(started_links.begin()+i);
}
i++;
}
}
}
paje.WriteEvents();
}
} // namespace ngcore
const char *header =
"%EventDef PajeDefineContainerType 0 \n"
"% Alias string \n"
"% Type string \n"
"% Name string \n"
"%EndEventDef \n"
"%EventDef PajeDefineVariableType 1 \n"
"% Alias string \n"
"% Type string \n"
"% Name string \n"
"% Color color \n"
"%EndEventDef \n"
"%EventDef PajeDefineStateType 2 \n"
"% Alias string \n"
"% Type string \n"
"% Name string \n"
"%EndEventDef \n"
"%EventDef PajeDefineEventType 3 \n"
"% Alias string \n"
"% Type string \n"
"% Name string \n"
"% Color color \n"
"%EndEventDef \n"
"%EventDef PajeDefineLinkType 4 \n"
"% Alias string \n"
"% Type string \n"
"% StartContainerType string \n"
"% EndContainerType string \n"
"% Name string \n"
"%EndEventDef \n"
"%EventDef PajeDefineEntityValue 5 \n"
"% Alias string \n"
"% Type string \n"
"% Name string \n"
"% Color color \n"
"%EndEventDef \n"
"%EventDef PajeCreateContainer 6 \n"
"% Time date \n"
"% Alias string \n"
"% Type string \n"
"% Container string \n"
"% Name string \n"
"%EndEventDef \n"
"%EventDef PajeDestroyContainer 7 \n"
"% Time date \n"
"% Type string \n"
"% Name string \n"
"%EndEventDef \n"
"%EventDef PajeSetVariable 8 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value double \n"
"%EndEventDef\n"
"%EventDef PajeAddVariable 9 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value double \n"
"%EndEventDef\n"
"%EventDef PajeSubVariable 10 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value double \n"
"%EndEventDef\n"
"%EventDef PajeSetState 11 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value string \n"
"%EndEventDef\n"
"%EventDef PajePushState 12 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value string \n"
"% Id string \n"
"%EndEventDef\n"
"%EventDef PajePopState 13 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"%EndEventDef\n"
"%EventDef PajeResetState 14 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"%EndEventDef\n"
"%EventDef PajeStartLink 15 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value string \n"
"% StartContainer string \n"
"% Key string \n"
"%EndEventDef\n"
"%EventDef PajeEndLink 16 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value string \n"
"% EndContainer string \n"
"% Key string \n"
"%EndEventDef\n"
"%EventDef PajeNewEvent 17 \n"
"% Time date \n"
"% Type string \n"
"% Container string \n"
"% Value string \n"
"%EndEventDef\n";

190
libsrc/core/paje_trace.hpp Normal file
View File

@ -0,0 +1,190 @@
#ifndef NETGEN_CORE_PAJE_TRACE_HPP
#define NETGEN_CORE_PAJE_TRACE_HPP
#include <limits>
#include <vector>
#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> 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<std::vector<Task> > tasks;
std::vector<Job> jobs;
std::vector<TimerEvent> timer_events;
std::vector<std::vector<ThreadLink> > 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

117
libsrc/core/profiler.cpp Normal file
View File

@ -0,0 +1,117 @@
#include <mutex>
#include "profiler.hpp"
namespace ngcore
{
std::vector<NgProfiler::TimerVal> 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<Logger> 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<std::mutex> 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

304
libsrc/core/profiler.hpp Normal file
View File

@ -0,0 +1,304 @@
#ifndef NETGEN_CORE_PROFILER_HPP
#define NETGEN_CORE_PROFILER_HPP
#include <chrono>
#include <string>
#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<TimerVal> timers;
NGCORE_API static TTimePoint * thread_times;
NGCORE_API static TTimePoint * thread_flops;
NGCORE_API static std::shared_ptr<Logger> 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<int>(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<typename TFunc>
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<double>::max();
int iteration = 0;
while(WallTime()<tend || iteration++ < min_iterations)
{
double t = -WallTime();
f();
t += WallTime();
tres = std::min(tres, t);
}
return tres;
}
} // namespace ngcore
#endif // NETGEN_CORE_PROFILER_HPP

30
libsrc/core/python_ngcore.cpp Normal file
View File

@ -0,0 +1,30 @@
#include <pybind11/pybind11.h>
#include "logging.hpp"
namespace py = pybind11;
using namespace ngcore;
PYBIND11_MODULE(pyngcore, m) // NOLINT
{
py::enum_<level::level_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.");
}

144
libsrc/core/symboltable.hpp Normal file
View File

@ -0,0 +1,144 @@
#ifndef NETGEN_CORE_SYMBOLTABLE_HPP
#define NETGEN_CORE_SYMBOLTABLE_HPP
#include <ostream>
#include <string>
#include <vector>
#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 T>
class SymbolTable
{
std::vector<std::string> names;
std::vector<T> data;
public:
using value_type = T;
using reference = typename std::vector<T>::reference;
using const_reference = typename std::vector<T>::const_reference;
/// Creates a symboltable
SymbolTable () = default;
SymbolTable (const SymbolTable<T> &) = default;
SymbolTable (SymbolTable<T> &&) noexcept = default;
~SymbolTable() = default;
SymbolTable& operator=(const SymbolTable<T>&) = default;
SymbolTable& operator=(SymbolTable<T>&&) = default;
template<typename T2=T>
auto DoArchive(Archive& ar) -> typename std::enable_if<is_archivable<T2>, 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<T>& Update(const SymbolTable<T>& tbl2)
{
for (size_t i = 0; i < tbl2.Size(); i++)
Set (tbl2.GetName(i), tbl2[i]);
return *this;
}
};
template <typename T>
std::ostream & operator<< (std::ostream & ost, const SymbolTable<T> & 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

42
libsrc/core/utils.cpp Normal file
View File

@ -0,0 +1,42 @@
#include "utils.hpp"
#include "logging.hpp"
#ifndef WIN32
#include <cxxabi.h>
#endif
#include <iostream>
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()<tend);
auto tick_end = GetTimeCounter();
tend = WallTime();
return (tick_end-tick_start)/(tend-tstart);
}();
const std::chrono::time_point<TClock> wall_time_start = TClock::now();
} // namespace ngcore

70
libsrc/core/utils.hpp Normal file
View File

@ -0,0 +1,70 @@
#ifndef NETGEN_CORE_UTILS_HPP
#define NETGEN_CORE_UTILS_HPP
#include <chrono>
#include <map>
#include <ostream>
#include <sstream>
#include <string>
#ifdef WIN32
#include <intrin.h> // for __rdtsc() CPU time step counter
#else
#include <x86intrin.h> // 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<TClock> wall_time_start;
// Time in seconds since program start
inline double WallTime () noexcept
{
std::chrono::time_point<TClock> now = TClock::now();
std::chrono::duration<double> 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 <class T>
inline std::string ToString (const T& t)
{
std::stringstream ss;
ss << t;
return ss.str();
}
template<typename T1, typename T2>
std::ostream& operator << (std::ostream& ost, const std::map<T1,T2>& map)
{
for(auto& val : map)
ost << "\n" << val.first << ": " << val.second;
return ost;
}
} // namespace ngcore
#endif // NETGEN_CORE_UTILS_HPP

6
libsrc/core/version.hpp
View File

@ -1,6 +1,7 @@
#ifndef NETGEN_CORE_VERSION_HPP
#define NETGEN_CORE_VERSION_HPP
#include <ostream>
#include <string>
#include <tuple>
@ -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

8
libsrc/csg/CMakeLists.txt
View File

@ -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 )

20
libsrc/csg/csgeom.cpp
View File

@ -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<int> surfind;
int i;

10
libsrc/csg/csgeom.hpp
View File

@ -102,7 +102,7 @@ namespace netgen
{
private:
/// all surfaces
SYMBOLTABLE<Surface*> surfaces;
SymbolTable<Surface*> surfaces;
public:
/// primitive of surface
@ -112,12 +112,12 @@ namespace netgen
Array<Surface*> delete_them;
/// all named solids
SYMBOLTABLE<Solid*> solids;
SymbolTable<Solid*> 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<TopLevelObject*> 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<Solid*> & GetSolids () const { return solids; }
const SymbolTable<Solid*> & GetSolids () const { return solids; }
void SetSplineCurve (const char * name, SplineGeometry<2> * spl);

11
libsrc/csg/python_csg.cpp
View File

@ -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<SPSolid> sol, shared_ptr<SPSolid> 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<SPSolid> s1,shared_ptr<SPSolid> s2, double factor)
{
auto singedge = new SingularEdge(1, -1, self, s1->GetSolid(), s2->GetSolid(), factor);

16
libsrc/csg/solid.cpp
View File

@ -419,9 +419,9 @@ namespace netgen
static Solid * CreateSolidExpr (istream & ist, const SYMBOLTABLE<Solid*> & solids);
static Solid * CreateSolidTerm (istream & ist, const SYMBOLTABLE<Solid*> & solids);
static Solid * CreateSolidPrim (istream & ist, const SYMBOLTABLE<Solid*> & solids);
static Solid * CreateSolidExpr (istream & ist, const SymbolTable<Solid*> & solids);
static Solid * CreateSolidTerm (istream & ist, const SymbolTable<Solid*> & solids);
static Solid * CreateSolidPrim (istream & ist, const SymbolTable<Solid*> & solids);
static void ReadString (istream & ist, char * str)
{
@ -461,7 +461,7 @@ namespace netgen
}
Solid * CreateSolidExpr (istream & ist, const SYMBOLTABLE<Solid*> & solids)
Solid * CreateSolidExpr (istream & ist, const SymbolTable<Solid*> & solids)
{
// cout << "create expr" << endl;
@ -484,7 +484,7 @@ namespace netgen
return s1;
}
Solid * CreateSolidTerm (istream & ist, const SYMBOLTABLE<Solid*> & solids)
Solid * CreateSolidTerm (istream & ist, const SymbolTable<Solid*> & solids)
{
// cout << "create term" << endl;
@ -508,7 +508,7 @@ namespace netgen
return s1;
}
Solid * CreateSolidPrim (istream & ist, const SYMBOLTABLE<Solid*> & solids)
Solid * CreateSolidPrim (istream & ist, const SymbolTable<Solid*> & 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<Solid*> & solids)
Solid * Solid :: CreateSolid (istream & ist, const SymbolTable<Solid*> & solids)
{
Solid * nsol = CreateSolidExpr (ist, solids);
nsol = new Solid (ROOT, nsol);

2
libsrc/csg/solid.hpp
View File

@ -158,7 +158,7 @@ namespace netgen
{ return maxh; }
void GetSolidData (ostream & ost, int first = 1) const;
static Solid * CreateSolid (istream & ist, const SYMBOLTABLE<Solid*> & solids);
static Solid * CreateSolid (istream & ist, const SymbolTable<Solid*> & solids);
static BlockAllocator ball;

22
libsrc/general/CMakeLists.txt
View File

@ -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
)

7
libsrc/general/array.hpp
View File

@ -243,10 +243,11 @@ namespace netgen
}
explicit Array(size_t asize)
: FlatArray<T, BASE, TIND> (asize, new T[asize])
: FlatArray<T, BASE, TIND> (asize, asize ? new T[asize] : nullptr)
{
allocsize = asize;
ownmem = 1;
allocsize = asize;
if(asize)
ownmem = 1;
}
/// Generate array in user data

20
libsrc/general/flags.cpp
View File

@ -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<double>&)numflags).Elem(name);
return & ((SymbolTable<double>&)numflags)[name];
else
return NULL;
}
@ -107,7 +107,7 @@ namespace netgen
double * Flags :: GetNumFlagPtr (const char * name)
{
if (numflags.Used (name))
return & ((SYMBOLTABLE<double>&)numflags).Elem(name);
return & ((SymbolTable<double>&)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<char*> 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<double> 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;
}

10
libsrc/general/flags.hpp
View File

@ -19,15 +19,15 @@ namespace netgen
class Flags
{
///
SYMBOLTABLE<char *> strflags;
SymbolTable<char *> strflags;
///
SYMBOLTABLE<double> numflags;
SymbolTable<double> numflags;
///
SYMBOLTABLE<int> defflags;
SymbolTable<int> defflags;
///
SYMBOLTABLE<Array<char*>*> strlistflags;
SymbolTable<Array<char*>*> strlistflags;
///
SYMBOLTABLE<Array<double>*> numlistflags;
SymbolTable<Array<double>*> numlistflags;
public:
///
DLL_HEADER Flags ();

168
libsrc/general/mpi_interface.hpp
View File

@ -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 <class T>
inline MPI_Datatype MyGetMPIType ( )
{ cerr << "ERROR in GetMPIType() -- no type found" << endl;return 0; }
template <>
inline MPI_Datatype MyGetMPIType<int> ( )
{ return MPI_INT; }
template <>
inline MPI_Datatype MyGetMPIType<double> ( )
{ return MPI_DOUBLE; }
template <>
inline MPI_Datatype MyGetMPIType<char> ( )
{ return MPI_CHAR; }
template<>
inline MPI_Datatype MyGetMPIType<size_t> ( )
{ return MPI_UINT64_T; }
#else
typedef int MPI_Datatype;
template <class T> inline MPI_Datatype MyGetMPIType ( ) { return 0; }
#endif
#ifdef PARALLEL
inline MPI_Comm MyMPI_SubCommunicator(MPI_Comm comm, Array<int> & 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<int> & 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 <class T>
MPI_Datatype MyGetMPIType ( )
{ cerr << "ERROR in GetMPIType() -- no type found" << endl;return 0; }
template <>
inline MPI_Datatype MyGetMPIType<int> ( )
{ return MPI_INT; }
template <>
inline MPI_Datatype MyGetMPIType<double> ( )
{ 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<char*> (s.c_str()), s.length(), MPI_CHAR, dest, tag, MPI_COMM_WORLD);
MPI_Send( const_cast<char*> (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 <class T, int BASE>
inline void MyMPI_Send (FlatArray<T, BASE> s, int dest, int tag)
inline void MyMPI_Send (FlatArray<T, BASE> s, int dest, int tag, MPI_Comm comm = ng_comm)
{
MPI_Send( &s.First(), s.Size(), MyGetMPIType<T>(), dest, tag, MPI_COMM_WORLD);
MPI_Send( &s.First(), s.Size(), MyGetMPIType<T>(), dest, tag, comm);
}
template <class T, int BASE>
inline void MyMPI_Recv ( FlatArray<T, BASE> s, int src, int tag)
inline void MyMPI_Recv ( FlatArray<T, BASE> s, int src, int tag, MPI_Comm comm = ng_comm)
{
MPI_Status status;
MPI_Recv( &s.First(), s.Size(), MyGetMPIType<T>(), src, tag, MPI_COMM_WORLD, &status);
MPI_Recv( &s.First(), s.Size(), MyGetMPIType<T>(), src, tag, comm, &status);
}
template <class T, int BASE>
inline void MyMPI_Recv ( Array <T, BASE> & s, int src, int tag)
inline void MyMPI_Recv ( Array <T, BASE> & 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<T>(), &len);
s.SetSize (len);
MPI_Recv( &s.First(), len, MyGetMPIType<T>(), src, tag, MPI_COMM_WORLD, &status);
MPI_Recv( &s.First(), len, MyGetMPIType<T>(), src, tag, comm, &status);
}
template <class T, int BASE>
inline int MyMPI_Recv ( Array <T, BASE> & s, int tag)
inline int MyMPI_Recv ( Array <T, BASE> & 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<T>(), &len);
s.SetSize (len);
MPI_Recv( &s.First(), len, MyGetMPIType<T>(), src, tag, MPI_COMM_WORLD, &status);
MPI_Recv( &s.First(), len, MyGetMPIType<T>(), src, tag, comm, &status);
return src;
}
@ -130,7 +201,7 @@ namespace netgen
*/
template <class T, int BASE>
inline MPI_Request MyMPI_ISend (FlatArray<T, BASE> s, int dest, int tag, MPI_Comm comm = MPI_COMM_WORLD)
inline MPI_Request MyMPI_ISend (FlatArray<T, BASE> s, int dest, int tag, MPI_Comm comm = ng_comm)
{
MPI_Request request;
MPI_Isend( &s.First(), s.Size(), MyGetMPIType<T>(), dest, tag, comm, &request);
@ -139,7 +210,7 @@ namespace netgen
template <class T, int BASE>
inline MPI_Request MyMPI_IRecv (FlatArray<T, BASE> s, int dest, int tag, MPI_Comm comm = MPI_COMM_WORLD)
inline MPI_Request MyMPI_IRecv (FlatArray<T, BASE> s, int dest, int tag, MPI_Comm comm = ng_comm)
{
MPI_Request request;
MPI_Irecv( &s.First(), s.Size(), MyGetMPIType<T>(), dest, tag, comm, &request);
@ -204,11 +275,10 @@ namespace netgen
template <typename T>
inline void MyMPI_ExchangeTable (TABLE<T> & send_data,
TABLE<T> & 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<int> send_sizes(ntasks);
Array<int> recv_sizes(ntasks);
@ -252,22 +322,22 @@ namespace netgen
template <class T>
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<T>(), 0, comm);
}
template <class T>
inline void MyMPI_Bcast (Array<T, 0> & s, MPI_Comm comm = MPI_COMM_WORLD)
inline void MyMPI_Bcast (Array<T, 0> & 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<T>(), 0, comm);
}
template <class T>
inline void MyMPI_Bcast (Array<T, 0> & s, int root, MPI_Comm comm = MPI_COMM_WORLD)
inline void MyMPI_Bcast (Array<T, 0> & s, int root, MPI_Comm comm = ng_comm)
{
int id;
MPI_Comm_rank(comm, &id);
@ -280,19 +350,19 @@ namespace netgen
}
template <class T, class T2>
inline void MyMPI_Allgather (const T & send, FlatArray<T2> recv, MPI_Comm comm)
inline void MyMPI_Allgather (const T & send, FlatArray<T2> recv, MPI_Comm comm = ng_comm)
{
MPI_Allgather( const_cast<T*> (&send), 1, MyGetMPIType<T>(), &recv[0], 1, MyGetMPIType<T2>(), comm);
}
template <class T, class T2>
inline void MyMPI_Alltoall (FlatArray<T> send, FlatArray<T2> recv, MPI_Comm comm)
inline void MyMPI_Alltoall (FlatArray<T> send, FlatArray<T2> recv, MPI_Comm comm = ng_comm)
{
MPI_Alltoall( &send[0], 1, MyGetMPIType<T>(), &recv[0], 1, MyGetMPIType<T2>(), comm);
}
// template <class T, class T2>
// inline void MyMPI_Alltoall_Block (FlatArray<T> send, FlatArray<T2> recv, int blocklen, MPI_Comm comm)
// inline void MyMPI_Alltoall_Block (FlatArray<T> send, FlatArray<T2> recv, int blocklen, MPI_Comm comm = ng_comm)
// {
// MPI_Alltoall( &send[0], blocklen, MyGetMPIType<T>(), &recv[0], blocklen, MyGetMPIType<T2>(), comm);
// }

4
libsrc/general/myadt.hpp
View File

@ -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"

5
libsrc/general/netgenout.hpp
View File

@ -4,14 +4,11 @@
// #include <ostream>
// #include <mystdlib.h>
// #include <meshing.hpp>
#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;

33
libsrc/general/ngexception.cpp
View File

@ -1,33 +0,0 @@
/**************************************************************************/
/* File: ngexception.cpp */
/* Author: Joachim Schoeberl */
/* Date: 16. Jan. 02 */
/**************************************************************************/
#include <myadt.hpp>
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;
}
}

34
libsrc/general/ngexception.hpp
View File

@ -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

12
libsrc/general/ngpython.hpp
View File

@ -3,6 +3,7 @@
#include <pybind11/pybind11.h>
#include <pybind11/operators.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>
namespace py = pybind11;
#include <iostream>
#include <sstream>
@ -66,16 +67,7 @@ namespace netgen
return static_cast<typename function_traits<Function>::pointer>(lambda);
}
template <class T>
inline std::string ToString (const T& t)
{
std::stringstream ss;
ss << t;
return ss.str();
}
}
} // namespace netgen
#endif

131
libsrc/general/profiler.cpp
View File

@ -1,131 +0,0 @@
/**************************************************************************/
/* File: profiler.cpp */
/* Author: Joachim Schoeberl */
/* Date: 19. Apr. 2002 */
/**************************************************************************/
#include <myadt.hpp>
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;
}

88
libsrc/general/profiler.hpp
View File

@ -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 <x86intrin.h> // 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<char*> (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<char*> (names[nr].c_str()));
}
#else
static void StartTimer (int nr)
{
starttimes[nr] = __rdtsc(); counts[nr]++;
// VT_USER_START (const_cast<char*> (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<char*> (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

52
libsrc/general/symbolta.cpp
View File

@ -1,52 +0,0 @@
/**************************************************************************/
/* File: symbolta.cc */
/* Author: Joachim Schoeberl */
/* Date: 01. Jun. 95 */
/**************************************************************************/
/*
Abstract data type Symbol Table
*/
#include <mystdlib.h>
#include <myadt.hpp>
#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

163
libsrc/general/symbolta.hpp
View File

@ -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 <char*> 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 T>
class SYMBOLTABLE : public BASE_SYMBOLTABLE
{
private:
/// Associated data
Array <T> 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<T> & operator= (SYMBOLTABLE<T> &);
};
template <class T>
inline SYMBOLTABLE<T> :: SYMBOLTABLE ()
{
;
}
template <class T>
inline INDEX SYMBOLTABLE<T> :: Size() const
{
return data.Size();
}
template <class T>
inline T & SYMBOLTABLE<T> :: Elem (const char * name)
{
int i = Index (name);
if (i)
return data.Elem (i);
else
return data.Elem(1);
}
template <class T>
inline const T & SYMBOLTABLE<T> :: Get (const char * name) const
{
int i;
i = Index (name);
if (i)
return data.Get(i);
else
return data.Get(1);
}
template <class T>
inline const T & SYMBOLTABLE<T> :: Get (int i) const
{
return data.Get(i);
}
template <class T>
inline const char* SYMBOLTABLE<T> :: GetName (int i) const
{
return names.Get(i);
}
template <class T>
inline void SYMBOLTABLE<T> :: 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 <class T>
inline bool SYMBOLTABLE<T> :: Used (const char * name) const
{
return (Index(name)) ? true : false;
}
template <class T>
inline void SYMBOLTABLE<T> :: DeleteAll ()
{
DelNames();
data.DeleteAll();
}
}
#endif

8
libsrc/geom2d/CMakeLists.txt
View File

@ -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)

7
libsrc/geom2d/genmesh2d.cpp
View File

@ -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++)

1
libsrc/geom2d/geometry2d.hpp
View File

@ -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);

49
libsrc/geom2d/python_geom2d.cpp
View File

@ -15,6 +15,22 @@ namespace netgen
DLL_HEADER void ExportGeom2d(py::module &m)
{
py::class_<SplineSegExt, shared_ptr<SplineSegExt>>
(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_<SplineGeometry2d, NetgenGeometry, shared_ptr<SplineGeometry2d>>
(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<SplineSegExt>(&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());

10
libsrc/gprim/CMakeLists.txt
View File

@ -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

6
libsrc/gprim/adtree.cpp
View File

@ -2400,13 +2400,13 @@ namespace netgen
Array<T> & 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);
}

6
libsrc/gprim/adtree.hpp
View File

@ -570,9 +570,9 @@ public:
{ tree->DeleteElement(pi); }
void GetIntersecting (const Point<dim> & pmin, const Point<dim> & pmax,
Array<T> & 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; };
};
}

3
libsrc/gprim/geomobjects.hpp
View File

@ -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 !
}
}

20
libsrc/include/nginterface.h
View File

@ -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);
///

91
libsrc/include/nginterface_v2.hpp
View File

@ -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 <typename T, int S>
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<Mesh> amesh = NULL);
void LoadMesh (const string & filename);
void LoadMesh (istream & str);
/** reuse a netgen-mesh **/
Ngx_Mesh (shared_ptr<Mesh> 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 <int DIM>
const Ng_Node<DIM> GetNode (int nr) const;
Ng_BufferMS<int,4> GetFaceEdges (int fnr) const;
template <int DIM>
int GetNNodes ();
@ -291,11 +336,17 @@ namespace netgen
// 3D only
// std::pair<int,int> GetBoundaryNeighbouringDomains (int bnr);
template <int DIM>
void SetRefinementFlag (size_t elnr, bool flag);
void Curve (int order);
void Refine (NG_REFINEMENT_TYPE reftype,
void (*taskmanager)(function<void(int,int)>) = &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<int,int*> GetDistantProcs (int nodetype, int locnum) const;
#endif
shared_ptr<Mesh> GetMesh () const { return mesh; }
shared_ptr<Mesh> 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;
};

7
libsrc/include/nginterface_v2_impl.hpp
View File

@ -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;
}

8
libsrc/interface/CMakeLists.txt
View File

@ -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

38
libsrc/interface/nginterface.cpp
View File

@ -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<char> 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);

211
libsrc/interface/nginterface_v2.cpp
View File

@ -31,39 +31,39 @@ namespace netgen
return hmesh;
}
Ngx_Mesh :: Ngx_Mesh (shared_ptr<Mesh> amesh)
{
if (amesh)
mesh = amesh;
else
mesh = netgen::mesh;
}
Ngx_Mesh :: Ngx_Mesh (shared_ptr<Mesh> 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<Mesh>();
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<Mesh>();
#ifdef PARALLEL
if (archive.Input()) {
mesh = make_shared<Mesh>();
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<int,4> Ngx_Mesh::GetFaceEdges (int fnr) const
{
const MeshTopology & topology = mesh->GetTopology();
ArrayMem<int,4> ia;
topology.GetFaceEdges (fnr+1, ia);
Ng_BufferMS<int,4> 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<void(int,int)>),
@ -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<Refinement&> (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<int,int*> Ngx_Mesh :: GetDistantProcs (int nodetype, int locnum) const
{
#ifdef PARALLEL
switch (nodetype)
{
case 0:
@ -1097,10 +1270,10 @@ namespace netgen
default:
return std::tuple<int,int*>(0,nullptr);
}
}
#else
return std::tuple<int,int*>(0,nullptr);
#endif
}
}

10
libsrc/linalg/CMakeLists.txt
View File

@ -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

16
libsrc/meshing/CMakeLists.txt
View File

@ -1,12 +1,4 @@
add_definitions(-DNGINTERFACE_EXPORTS)
if(NOT WIN32)
set(mesh_object_libs
$<TARGET_OBJECTS:la>
$<TARGET_OBJECTS:gprim>
$<TARGET_OBJECTS:gen>
)
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

1
libsrc/meshing/adfront3.hpp
View File

@ -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;
///

3
libsrc/meshing/bisect.cpp
View File

@ -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++)

273
libsrc/meshing/curvedelems.cpp
View File

@ -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<int> 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<double>)) 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<double>)) 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;

316
libsrc/meshing/delaunay.cpp
View File

@ -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<int> & 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<Element> newels;
// Array<Element> newels;
Array<DelaunayTet> 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<DelaunayTet> & tempels,
int oldnp, DelaunayTet & startel, Point3d & pmin, Point3d & pmax)
{
Array<Point<3> > centers;
static Timer t("Meshing3::Delaunay1"); RegionTimer reg(t);
static Timer tloop("Meshing3::Delaunay1 loop");
Array<Point<3>> centers;
Array<double> 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<PointIndex::BASE> 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<int> freelist;
int cntp = 0;
MeshNB meshnb (tempels, mesh.GetNP() + 5);
@ -689,13 +595,12 @@ namespace netgen
list.AddElement (1);
Array<int> 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<int> 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<int,PointIndex::BASE> 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<INDEX_2> 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))

7
libsrc/meshing/global.cpp
View File

@ -31,6 +31,9 @@ namespace netgen
DLL_HEADER shared_ptr<NetgenGeometry> ng_geometry;
// TraceGlobal glob2("global2");
// global communicator for netgen
DLL_HEADER MPI_Comm ng_comm = MPI_COMM_WORLD;
weak_ptr<Mesh> global_mesh;
void SetGlobalMesh (shared_ptr<Mesh> 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)

4
libsrc/meshing/global.hpp
View File

@ -59,6 +59,10 @@ namespace netgen
DLL_HEADER extern weak_ptr<Mesh> global_mesh;
DLL_HEADER void SetGlobalMesh (shared_ptr<Mesh> m);
// global communicator for netgen (dummy if no MPI)
extern DLL_HEADER MPI_Comm ng_comm;
}
#endif

78
libsrc/meshing/hprefinement.cpp
View File

@ -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; i<np ; i++)
pnums[i] = el[i];
index = el.GetIndex();
const Point3d * points =
MeshTopology :: GetVertices (el.GetType());
for(int i=0;i<np;i++)
for(int l=0;l<3;l++)
param[i][l] = points[i].X(l+1);
type = HP_NONE;
domin=-1; domout=-1; // he: needed for segments
}
HPRefElement :: HPRefElement(Element2d & el)
HPRefElement :: HPRefElement(Element2d & el) :
levelx(0), levely(0), levelz(0), type(HP_NONE), index(el.GetIndex()), np(el.GetNV()), domin(-1), domout(-1) //domin,out for segements
{
//Reset();
np = el.GetNV();
for (int i=0; i<np ; i++)
pnums[i] = el[i];
index = el.GetIndex();
const Point3d * points =
MeshTopology :: GetVertices (el.GetType());
for(int i=0;i<np;i++)
for(int l=0;l<3;l++)
param[i][l] = points[i].X(l+1);
type = HP_NONE;
domin=-1; domout=-1; // he: needed for segments
}
HPRefElement :: HPRefElement(Segment & el)
HPRefElement :: HPRefElement(Segment & el) :
levelx(0), levely(0), levelz(0), type(HP_NONE), np(2), domin(el.domin), domout(el.domout), singedge_left(el.singedge_left), singedge_right(el.singedge_right)
{
//Reset();
np = 2;
for (int i=0; i<np ; i++)
pnums[i] = el[i];
const Point3d * points =
@ -86,35 +81,18 @@ namespace netgen
param[i][1] = -1; param[i][2] = -1;
}
*/
singedge_left = el.singedge_left;
singedge_right = el.singedge_right;
type = HP_NONE;
// he: needed for orientation!
domin = el.domin;
domout = el.domout;
}
HPRefElement :: HPRefElement(HPRefElement & el)
HPRefElement :: HPRefElement(HPRefElement & el) :
np(el.np), levelx(el.levelx), levely(el.levely), levelz(el.levelz), type(el.type), domin(el.domin), domout(el.domout), index(el.index), coarse_elnr(el.coarse_elnr), singedge_left(el.singedge_left), singedge_right(el.singedge_right)
{
//Reset();
np = el.np;
for (int i=0; i<np ; i++)
{
pnums[i] = el[i];
for(int l=0; l<3; l++) param[i][l] = el.param[i][l];
}
index = el.index;
levelx = el.levelx;
levely = el.levely;
levelz = el.levelz;
type = el.type;
coarse_elnr = el.coarse_elnr;
singedge_left = el.singedge_left;
singedge_right = el.singedge_right;
domin = el.domin; // he: needed for segments
domout=el.domout;
}
}
void HPRefElement :: SetType (HPREF_ELEMENT_TYPE t)
@ -706,8 +684,7 @@ namespace netgen
{
HPRefElement el = elements[i];
HPRef_Struct * hprs = Get_HPRef_Struct (el.type);
int newlevel = el.levelx + 1;
int newlevel = el.levelx;
int oldnp = 0;
switch (hprs->geom)
{
@ -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<HPRefElement> & hpelements = *mesh.hpelements;
InitHPElements(mesh,hpelements);
Array<int> nplevel;
nplevel.Append (mesh.GetNP());

8
libsrc/meshing/improve3.cpp
View File

@ -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);

4
libsrc/meshing/meshclass.cpp
View File

@ -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;

26
libsrc/meshing/meshclass.hpp
View File

@ -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<Segment, 0, size_t> segments;
@ -132,8 +130,8 @@ namespace netgen
/// mesh access semaphors.
NgMutex majormutex;
SYMBOLTABLE< Array<int>* > userdata_int;
SYMBOLTABLE< Array<double>* > userdata_double;
SymbolTable< Array<int>* > userdata_int;
SymbolTable< Array<double>* > 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<int> & neloc) const;
/// loads a mesh sent from master processor
void ReceiveParallelMesh ();
#else
void Distribute () {}
void SendRecvMesh () {}
void Distribute (Array<int> & volume_weights, Array<int> & surface_weights,
Array<int> & segment_weights){ }
#endif

6
libsrc/meshing/meshfunc.cpp
View File

@ -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;

49
libsrc/meshing/meshing3.cpp
View File

@ -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<Point3d, PointIndex::BASE> 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;

92
libsrc/meshing/meshtype.cpp
View File

@ -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");
}

25
libsrc/meshing/meshtype.hpp
View File

@ -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)

176
libsrc/meshing/parallelmesh.cpp
View File

@ -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<MPI_Request> 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<MeshTopology&>(GetTopology()).Update();
PrintMessage ( 3, "Sending nr of elements");
MPI_Comm comm = this->comm;
Array<int> 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<PointIndex> 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<MPI_Request> 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<osegs_both.Size(); l++) {
int segnp = (*this)[osegs_both[l]].GetNP();
if(segnp!=segnp)
int segnp2 = (*this)[osegs_both[l]].GetNP();
if(segnp!=segnp2)
throw NgException("Tried to identify non-curved and curved Segment!");
}
for(int l = 0; l<osegs_both.Size(); l++) {
@ -676,7 +677,7 @@ namespace netgen
});
// distrubute segment data
for (int dest = 1; dest < ntasks; dest++)
sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5));
sendrequests.Append (MyMPI_ISend(segm_buf[dest], dest, MPI_TAG_MESH+5, comm));
PrintMessage ( 3, "now wait ...");
@ -700,9 +701,9 @@ namespace netgen
compiled_bcnames[tot_bcsize++] = (*bcnames[k])[j];
for(int k=1;k<ntasks;k++) {
sendrequests[6*(k-1)] = MyMPI_ISend(FlatArray<int>(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<int>(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<mat_sizes[k];j++)
compiled_mats[tot_matsize++] = (*materials[k])[j];
for(int k=1;k<ntasks;k++) {
sendrequests[6*(k-1)+3] = MyMPI_ISend(FlatArray<int>(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<int>(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<int> 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<int> 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<int> 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<int> 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<double> 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<int> 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<double> 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<int> 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<nbcs;k++)
size_bc += bcs[k];
char compiled_bcnames[size_bc];
MPI_Recv(compiled_bcnames, size_bc, MPI_CHAR, 0, MPI_TAG_MESH+6, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(compiled_bcnames, size_bc, MPI_CHAR, 0, MPI_TAG_MESH+6, comm, MPI_STATUS_IGNORE);
SetNBCNames(nbcs);
int cnt = 0;
@ -957,14 +959,14 @@ namespace netgen
/** Recv mat-names **/
int nmats;
MyMPI_Recv(nmats, 0, MPI_TAG_MESH+6);
MyMPI_Recv(nmats, 0, MPI_TAG_MESH+6, comm);
Array<int> 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<nmats;k++)
size_mats += matsz[k];
char compiled_mats[size_mats];
MPI_Recv(compiled_mats, size_mats, MPI_CHAR, 0, MPI_TAG_MESH+6, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(compiled_mats, size_mats, MPI_CHAR, 0, MPI_TAG_MESH+6, comm, MPI_STATUS_IGNORE);
cnt = 0;
materials.SetSize(nmats);
for(int k=0;k<nmats;k++) {
@ -973,7 +975,7 @@ namespace netgen
cnt += matsz[k];
}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Barrier(comm);
int timerloc = NgProfiler::CreateTimer ("Update local mesh");
int timerloc2 = NgProfiler::CreateTimer ("CalcSurfacesOfNode");
@ -982,7 +984,8 @@ namespace netgen
stringstream str;
str << "p" << id << ": got " << GetNE() << " elements and "
<< GetNSE() << " surface elements";
cout << str.str() << endl;
PrintMessage(2, str.str());
// cout << str.str() << endl;
// PrintMessage (2, "Got ", GetNE(), " elements and ", GetNSE(), " surface elements");
// PrintMessage (2, "Got ", GetNSE(), " surface elements");
@ -1008,8 +1011,9 @@ namespace netgen
// call it only for the master !
void Mesh :: Distribute ()
{
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;
@ -1068,30 +1072,46 @@ namespace netgen
eptr.Append (eind.Size());
Array<idx_t> 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<bool, PointIndex::BASE> 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<int, PointIndex::BASE> pnt2el(cnt);
loop_els([&](auto vertex, int index)
{
@ -1277,8 +1293,9 @@ namespace netgen
// call it only for the master !
void Mesh :: Distribute (Array<int> & volume_weights , Array<int> & surface_weights, Array<int> & 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<idx_t> 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;

15
libsrc/meshing/paralleltop.cpp
View File

@ -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;

165
libsrc/meshing/python_mesh.cpp
View File

@ -17,6 +17,22 @@ namespace netgen
{
extern bool netgen_executable_started;
extern shared_ptr<NetgenGeometry> ng_geometry;
#ifdef PARALLEL
/** we need allreduce in python-wrapped communicators **/
template <typename T>
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<T>(), op, comm);
return global_d;
}
#else
enum { MPI_SUM = 0, MPI_MIN = 1, MPI_MAX = 2 };
typedef int MPI_Op;
template <typename T>
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_<Point<3>> (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_<Vec<3>> (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_<Element>(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<PointIndex>(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<PointIndex>(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<PointIndex>(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<PointIndex>(vertices[i])();
newel->SetIndex(index);
}
else
throw NgException ("cannot create element");
std::map<int, ELEMENT_TYPE> 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<np; i++)
(*newel)[i] = py::cast<PointIndex>(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<PointIndex>(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<PointIndex>(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_<Mesh,shared_ptr<Mesh>>(m, "Mesh")
// .def(py::init<>("create empty mesh"))
.def(py::init( [] (int dim)
.def(py::init( [] (int dim, shared_ptr<PyMPI_Comm> pycomm)
{
auto mesh = make_shared<Mesh>();
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<Mesh>())
.def_property_readonly("comm", [](const Mesh & amesh)
{ return make_shared<PyMPI_Comm>(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<Mesh>)
.def_property_readonly("_timestamp", &Mesh::GetTimeStamp)
.def("Distribute", [](shared_ptr<Mesh> self, shared_ptr<PyMPI_Comm> 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<PyMPI_Comm> pycomm) {
auto mesh = make_shared<Mesh>();
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_<PyMPI_Comm, shared_ptr<PyMPI_Comm>> (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<int> proc_list) {
Array<int> 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<PyMPI_Comm>(subcomm, true);
/*
Array<int> procs;
if (py::extract<py::list> (proc_list).check()) {
py::list pylist = py::extract<py::list> (proc_list)();
procs.SetSize(py::len(pyplist));
for (int i = 0; i < py::len(pylist); i++)
procs[i] = py::extract<int>(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<PyMPI_Comm>(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<PyMPI_Comm>(subcomm, true);
*/
}, py::arg("procs"));
;
}
PYBIND11_MODULE(libmesh, m) {

62
libsrc/meshing/secondorder.cpp
View File

@ -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()));
}

2
libsrc/meshing/smoothing3.cpp
View File

@ -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";

16
libsrc/meshing/topology.cpp
View File

@ -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<int> 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<int> 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<int> 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<int> vert2face(2*max_face_on_vertex+10);

40
libsrc/meshing/topology.hpp
View File

@ -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:

6
libsrc/occ/CMakeLists.txt
View File

@ -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)

9
libsrc/occ/occmeshsurf.cpp
View File

@ -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);
}

6
libsrc/stlgeom/CMakeLists.txt
View File

@ -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

1
libsrc/stlgeom/meshstlsurface.cpp
View File

@ -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;

6
libsrc/visualization/CMakeLists.txt
View File

@ -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

10
libsrc/visualization/soldata.hpp
View File

@ -1,7 +1,7 @@
#ifndef FILE_SOLDATA
#define FILE_SOLDATA
#include <myadt.hpp> // 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<double>" << endl;
return false;

2
libsrc/visualization/vsmesh.cpp
View File

@ -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;

15
libsrc/visualization/vssolution.cpp
View File

@ -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);

6
ng/ngpkg.cpp
View File

@ -1840,9 +1840,9 @@ namespace netgen
SYMBOLTABLE<VisualScene*> & GetVisualizationScenes ()
SymbolTable<VisualScene*> & GetVisualizationScenes ()
{
static SYMBOLTABLE<VisualScene*> vss;
static SymbolTable<VisualScene*> 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)
{

4
nglib/CMakeLists.txt
View File

@ -8,9 +8,6 @@ if(WIN32)
$<TARGET_OBJECTS:geom2d>
$<TARGET_OBJECTS:csg>
$<TARGET_OBJECTS:stl>
$<TARGET_OBJECTS:gen>
$<TARGET_OBJECTS:la>
$<TARGET_OBJECTS:gprim>
$<TARGET_OBJECTS:visual>
$<TARGET_OBJECTS:occ>
@ -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} )

4
nglib/nglib.cpp
View File

@ -1211,14 +1211,14 @@ namespace netgen
}
*/
/*
#ifndef WIN32
void ResetTime ()
{
;
}
#endif
*/
void MyBeep (int i)

98
python/geom2d.py
View File

@ -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']

108
python/read_gmsh.py
View File

@ -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

8
tests/build_nospdlog.sh Normal file
View File

@ -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

19
tests/catch/CMakeLists.txt
View File

@ -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)

38
tests/catch/archive.cpp
View File

@ -236,6 +236,34 @@ void testMultipleInheritance(Archive& in, Archive& out)
}
}
void testMap(Archive& in, Archive& out)
{
map<string, VersionInfo> map1;
map1["netgen"] = "v6.2.1901";
out & map1;
out.FlushBuffer();
map<string, VersionInfo> 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<stringstream>();
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<stringstream>();
TextOutArchive out(stream);
TextInArchive in(stream);

64
tests/catch/symboltable.cpp Normal file
View File

@ -0,0 +1,64 @@
#include "catch.hpp"
#include <../core/ngcore.hpp>
using namespace ngcore;
using namespace std;
TEST_CASE("Symboltable")
{
SymbolTable<int> 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<int> 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<std::stringstream>();
BinaryOutArchive ao(ss);
ao & table;
ao.FlushBuffer();
BinaryInArchive ai(ss);
SymbolTable<int> read;
ai & read;
for(size_t i = 0; i<table.Size(); i++)
{
CHECK(read[i] == table[i]);
CHECK(read.GetName(i) == table.GetName(i));
}
table.DeleteAll();
CHECK(table.Size() == 0);
// SymbolTable<bool> is special because of vector<bool> is special...
SymbolTable<bool> 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<bool> bread;
ai & bread;
CHECK(bread["true"]);
CHECK(!bread["false"]);
}