diff --git a/CMake/FindTBB.cmake b/CMake/FindTBB.cmake
index 71565cc7a05e3165877c0234a4d4cd2e038693a2..6f17524cea41116149174a1995750f01e1254ffc 100644
--- a/CMake/FindTBB.cmake
+++ b/CMake/FindTBB.cmake
@@ -78,7 +78,7 @@ SET(TBB_CXX_FLAGS_DEBUG "-DTBB_USE_THREADING_TOOLS=1 -DTBB_USE_DEBUG=1")
# --- List of libraries to be found ---
-IF(DEFINED ENV{ONEAPI_ROOT})
+IF(DEFINED ENV{ONEAPI_ROOT} AND MSVC)
set(_TBB_TBBLIB_NAME "tbb12")
ELSE()
set(_TBB_TBBLIB_NAME "tbb")
@@ -245,37 +245,37 @@ unset(_TBB_INCLUDE_DIRS CACHE)
# Versions
# ----------------------------------------------------------------------------
-SET(_TBB_VERSION_INTERFACE ${TBB_VERSION_INTERFACE})
-IF( (NOT TBB_VERSION_INTERFACE) AND TBB_INCLUDE_DIRS AND NOT DEFINED ENV{ONEAPI_ROOT})
- FILE(READ "${TBB_INCLUDE_DIRS}/tbb/tbb_stddef.h" _FILE)
- set(_TBB_VERSION_MAJOR 0)
- set(_TBB_VERSION_MINOR 0)
- STRING(REGEX REPLACE ".*#define TBB_INTERFACE_VERSION ([0-9]+).*" "\\1" _TBB_VERSION_INTERFACE "${_FILE}")
+# Find file containing version info and read the version
+IF(EXISTS "${TBB_INCLUDE_DIRS}/tbb/tbb_stddef.h")
+ FILE(READ "${TBB_INCLUDE_DIRS}/tbb/tbb_stddef.h" _FILE) # older TBB
STRING(REGEX REPLACE ".*#define TBB_VERSION_MAJOR ([0-9]+).*" "\\1" _TBB_VERSION_MAJOR "${_FILE}")
STRING(REGEX REPLACE ".*#define TBB_VERSION_MINOR ([0-9]+).*" "\\1" _TBB_VERSION_MINOR "${_FILE}")
- set(_TBB_VERSION_STRING "${_TBB_VERSION_MAJOR}.${_TBB_VERSION_MINOR}")
- unset(_FILE)
+ SET(TBB_VERSION ${_TBB_VERSION_MAJOR}.${_TBB_VERSION_MINOR})
+ELSEIF(EXISTS "${TBB_INCLUDE_DIRS}/oneapi/tbb/version.h")
+ FILE(READ "${TBB_INCLUDE_DIRS}/oneapi/tbb/version.h" _FILE) # newer oneTBB
+ STRING(REGEX REPLACE ".*#define TBB_VERSION_MAJOR ([0-9]+).*" "\\1" _TBB_VERSION_MAJOR "${_FILE}")
+ STRING(REGEX REPLACE ".*#define TBB_VERSION_MINOR ([0-9]+).*" "\\1" _TBB_VERSION_MINOR "${_FILE}")
+ STRING(REGEX REPLACE ".*#define TBB_VERSION_PATCH ([0-9]+).*" "\\1" _TBB_VERSION_PATCH "${_FILE}")
+ SET(TBB_VERSION ${_TBB_VERSION_MAJOR}.${_TBB_VERSION_MINOR}.${_TBB_VERSION_PATCH})
+ELSE()
+ MESSAGE(STATUS "TBB_LIBRARIES: ${TBB_LIBRARIES}")
+ MESSAGE(STATUS "TBB_INCLUDE_DIRS: ${TBB_INCLUDE_DIRS}")
+ MESSAGE(FATAL_ERROR "Cannot find version file in ${TBB_INCLUDE_DIRS}/tbb.")
ENDIF()
-IF(NOT TBB_VERSION_INTERFACE)
- SET(TBB_VERSION_INTERFACE ${_TBB_VERSION_INTERFACE} CACHE STRING "" FORCE)
- SET(TBB_VERSION_MAJOR ${_TBB_VERSION_MAJOR} CACHE INTERNAL "" FORCE)
- SET(TBB_VERSION_MINOR ${_TBB_VERSION_MINOR} CACHE INTERNAL "" FORCE)
- SET(TBB_VERSION_STRING ${_TBB_VERSION_STRING} CACHE STRING "" FORCE)
- mark_as_advanced(TBB_VERSION_INTERFACE)
- mark_as_advanced(TBB_VERSION_MAJOR)
- mark_as_advanced(TBB_VERSION_MINOR)
- mark_as_advanced(TBB_VERSION_STRING)
+# If file has been found, also read the interface version
+STRING(REGEX REPLACE ".*#define TBB_INTERFACE_VERSION ([0-9]+).*" "\\1" _TBB_INTERFACE_VERSION "${_FILE}")
+SET(TBB_INTERFACE_VERSION ${_TBB_INTERFACE_VERSION} CACHE STRING "" FORCE)
+IF (_VERB)
+ MESSAGE(STATUS "TBB version: ${TBB_VERSION} (${TBB_INTERFACE_VERSION})")
ENDIF()
# ----------------------------------------------------------------------------
include(FindPackageHandleStandardArgs)
# handle the QUIETLY and REQUIRED arguments and set TBB_FOUND to TRUE
# if all listed variables are TRUE
-find_package_handle_standard_args(TBB DEFAULT_MSG
- TBB_LIBRARIES TBB_INCLUDE_DIRS
- TBB_LIB_PATH TBB_LIB_PATH_DEBUG)
-
-#MESSAGE(STATUS "TBB_FOUND = ${TBB_FOUND}")
+find_package_handle_standard_args(TBB
+ REQUIRED_VARS TBB_LIBRARIES TBB_INCLUDE_DIRS TBB_LIB_PATH TBB_LIB_PATH_DEBUG
+ VERSION_VAR TBB_VERSION)
# ----------------------------------------------------------------------------
unset(_VERB)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5b5da42d071ae94d6989d1bb93a7538f3467114a..8bb8eca85ec84335a03f532aede76d2d355420f3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,11 +15,7 @@
# ----------------------------------------------------------------------------
PROJECT(amfe)
# ----------------------------------------------------------------------------
-CMAKE_MINIMUM_REQUIRED(VERSION 3.1)
-IF(${CMAKE_VERSION} VERSION_GREATER "3.14.0") # we might want to update the project and use the NEW behavior here
- cmake_policy(SET CMP0078 OLD)
- cmake_policy(SET CMP0086 OLD)
-ENDIF()
+CMAKE_MINIMUM_REQUIRED(VERSION 3.14)
# -- I/O
# Lib/Exe dir
@@ -66,17 +62,12 @@ OPTION(USE_MUMPS "Compile with MUMPS" ON)
# -- DEPENDENCIES
# Python
-IF (CMAKE_VERSION VERSION_LESS 3.12.0)
- FIND_PACKAGE(PythonInterp 3.6 REQUIRED)
- FIND_PACKAGE(PythonLibs 3.6 REQUIRED)
-ELSE()
- FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
- # use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
- SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
- SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
- SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS})
- SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})
-ENDIF()
+# use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
+FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
+SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
+SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
+SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS})
+SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})
MESSAGE(STATUS "PYTHON_EXECUTABLE:FILEPATH=${PYTHON_EXECUTABLE}")
MESSAGE(STATUS "PYTHON_LIBRARIES:FILEPATH=${PYTHON_LIBRARIES}")
MESSAGE(STATUS "PYTHON_INCLUDE_PATH:FILEPATH=${PYTHON_INCLUDE_PATH}")
diff --git a/fwk/_src/CMakeLists.txt b/fwk/_src/CMakeLists.txt
index 78c8658478fbd983c3d1625165f0e6957d27a216..fa26e7c3cce7349259a4e954a2738d4f5fccb86a 100644
--- a/fwk/_src/CMakeLists.txt
+++ b/fwk/_src/CMakeLists.txt
@@ -19,29 +19,16 @@ INCLUDE(${SWIG_USE_FILE})
FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
FILE(GLOB ISRCS *.i)
-SET(CMAKE_SWIG_FLAGS "")
SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} "-interface" "_fwkw") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(fwkw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET fwkw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+MACRO_DebugPostfix(fwkw)
-SET(SWINCFLAGS
--I${PROJECT_SOURCE_DIR}/fwk/src
-)
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
-
-if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
- SWIG_ADD_MODULE(fwkw python ${ISRCS} ${SRCS})
-else()
- SWIG_ADD_LIBRARY(fwkw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-endif()
-MACRO_DebugPostfix(_fwkw)
-
-TARGET_INCLUDE_DIRECTORIES(_fwkw PRIVATE ${PROJECT_SOURCE_DIR}/fwk/_src
- ${PYTHON_INCLUDE_PATH}
-)
-
-SWIG_LINK_LIBRARIES(fwkw
- fwk ${PYTHON_LIBRARIES}
-)
+TARGET_INCLUDE_DIRECTORIES(fwkw PRIVATE ${PROJECT_SOURCE_DIR}/fwk/_src
+ ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(fwkw PRIVATE fwk ${PYTHON_LIBRARIES})
INSTALL(FILES "${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/fwkw.py"
DESTINATION ${CMAKE_INSTALL_PREFIX})
-INSTALL(TARGETS _fwkw DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(TARGETS fwkw DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/fwk/wutils.py b/fwk/wutils.py
index d7e2a7a0a6c80da61aa796c66e5f787bed6bd36d..d5baa08e243c06eb43a2344f9314a93f4463408c 100644
--- a/fwk/wutils.py
+++ b/fwk/wutils.py
@@ -1,13 +1,13 @@
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
-#
+#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
-#
+#
# http://www.apache.org/licenses/LICENSE-2.0
-#
+#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -96,7 +96,7 @@ def findbins(modname, dirB=['build','wavesB'], verb=False): # changer wavesB en
print("[findbins] looking for linux/mingw libs")
libpath.append( os.path.join(spath, 'bin') ) # also for win-mingw
- # check the paths one by one
+ # check the paths one by one
for p in libpath:
if os.path.isdir(p):
if verb:
@@ -133,40 +133,35 @@ def initDLL():
# ------------------------------------------------------------------------------
-def initMKL(nthreads, verb=False):
+def initMKL(verb=False):
+ """Initialize the environment for MKL
+ Since we use the SDL for MKL, we need to configure it at runtime.
+ This is MUST be performed before any calls to MKL.
+ Could be improved by performing the init from the C++ and actually checking that the threading layer is the one we want, instead of using env var.
+ """
import os
- # we try to have full control over the threading environment
- # => set MKL/OMP variables if (and only if) they are not set!
- # number of "default" threads for MKL (available MKL domains: ALL, BLAS, VML, PARDISO)
- if not 'MKL_NUM_THREADS' in os.environ:
- os.environ['MKL_NUM_THREADS'] = '%d' % nthreads
- os.environ['MKL_DOMAIN_NUM_THREADS'] = 'MKL_DOMAIN_ALL=%d' % nthreads
- # threading interface
- if not 'MKL_THREADING_LAYER' in os.environ:
- if nthreads == 1:
- os.environ['MKL_THREADING_LAYER'] = 'SEQUENTIAL'
- else:
- os.environ['MKL_THREADING_LAYER'] = 'TBB'
- # interface (32 or 64 bits)
+ # Define threading interface
if not 'MKL_INTERFACE_LAYER' in os.environ:
os.environ['MKL_INTERFACE_LAYER'] = 'LP64' # 32 bits (default), ILP64 for 64 bits
- # we do not want to use OpenMP
- if not 'OMP_NUM_THREADS' in os.environ:
- os.environ['OMP_NUM_THREADS'] = '1'
+ # Define threading layer
+ if not 'MKL_THREADING_LAYER' in os.environ:
+ os.environ['MKL_THREADING_LAYER'] = 'TBB' # sequential computation will be performed with 1 TBB thread
+ # Force number of OMP threads to 1, TBB threads are controlled from inside the C++ code
+ if not 'MKL_NUM_THREADS' in os.environ:
+ os.environ['MKL_NUM_THREADS'] = '1'
+ os.environ['MKL_DOMAIN_NUM_THREADS'] = 'MKL_DOMAIN_ALL=1'
+ # Force number of OMP threads to remain unchanged
if not 'MKL_DYNAMIC' in os.environ:
os.environ['MKL_DYNAMIC'] = 'FALSE'
- if not 'OMP_DYNAMIC' in os.environ:
- os.environ['OMP_DYNAMIC'] = 'FALSE'
- # display variables
- print(f'Envrionment initialized for MKL with {nthreads} threads.')
+ # Display variables
if verb:
- for s in ['OMP_', 'MKL', 'KMP']:
- print('* %s environment:' % s)
+ for s in ['OMP', 'MKL', 'KMP']:
+ print('* %s environment:' % s)
for key, value in os.environ.items():
if s in key:
print(' %s = %s' % (key, value))
print('')
-
+
# ------------------------------------------------------------------------------
def setupwdir(testname):
@@ -188,7 +183,7 @@ def setupwdir(testname):
name = "noname"
import os, os.path
- # builds the name of the workspace folder
+ # builds the name of the workspace folder
#print "__file__=",__file__
dir1 = os.path.abspath(os.path.dirname(__file__)+os.sep+"..")+os.sep
print("dir1=",dir1)
@@ -198,7 +193,7 @@ def setupwdir(testname):
resdir = testname[len(common):].replace(os.sep,"_")
resdir = os.path.splitext(resdir)[0] # remove ".py"
if siz>1:
- # append the mpi size if mpi is used
+ # append the mpi size if mpi is used
# (so that the same test can be run concurrently with
# sevral mpi sizes in ctest)
resdir += "_mpi%d" % siz
@@ -220,10 +215,10 @@ def setupwdir(testname):
os.remove(fpth)
elif os.path.isdir(fpth):
shutil.rmtree(fpth)
-
+
if siz>1: # if multiple process => send sync to slaves.
comm.Barrier()
-
+
# we can now safely go into the new folder
os.chdir(wdir)
@@ -283,7 +278,7 @@ def run():
# parse args
args = parseargs()
- initMKL(args.k, verb=True) # initialize threading layer and number of threads
+ initMKL(verb=True) # initialize MKL environment (threading interface, layer and number of threads)
if args.fpe:
fwk.enableFpe() # enables floating point exceptions at runtime
print(ccolors.ANSI_YELLOW + 'Floating point exceptions will cause numpy to overflow!' + ccolors.ANSI_RESET)
@@ -305,4 +300,5 @@ def run():
print("starting test", testname)
print("time:", time.strftime("%c"))
print("hostname:", socket.gethostname())
- exec(open(testname, 'r', encoding='utf8').read(), globals(), globals()) # exec at global level!!
+ script = open(testname, 'r', encoding='utf-8').read()
+ exec(compile(script, testname, 'exec'), {'__file__': testname, '__name__':'__main__'}) # compile to get clean backtrace
diff --git a/init.py b/init.py
index d6bb154c89a0d4bc6e06f5fb34ca775c16b8e303..1d1037f1a20bdaa95e7094b807f50f1d98fed800 100644
--- a/init.py
+++ b/init.py
@@ -2,3 +2,6 @@
# Add this dir to the python path so that sub-dir can be imported as sub-module
import os.path, sys
sys.path.append(os.path.dirname(os.path.realpath(__file__)))
+# Initialize MKL environment
+from .fwk import wutils
+wutils.initMKL()
diff --git a/tbox/_src/CMakeLists.txt b/tbox/_src/CMakeLists.txt
index 1dea57d843faf64c93ee41ba2a9ce4ceee12ea56..9e309715d7fd92b4cd209b1afd2e37159a8b2435 100644
--- a/tbox/_src/CMakeLists.txt
+++ b/tbox/_src/CMakeLists.txt
@@ -19,33 +19,17 @@ INCLUDE(${SWIG_USE_FILE})
FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
FILE(GLOB ISRCS *.i)
-SET(CMAKE_SWIG_FLAGS "")
SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} "-interface" "_tboxw") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(tboxw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET tboxw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+MACRO_DebugPostfix(tboxw)
-SET(SWINCFLAGS
--I${PROJECT_SOURCE_DIR}/tbox/src
--I${PROJECT_SOURCE_DIR}/tbox/_src
--I${PROJECT_SOURCE_DIR}/fwk/src
--I${PROJECT_SOURCE_DIR}/fwk/_src
-)
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
-
-if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
- SWIG_ADD_MODULE(tboxw python ${ISRCS} ${SRCS})
-else()
- SWIG_ADD_LIBRARY(tboxw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-endif()
-MACRO_DebugPostfix(_tboxw)
-
-TARGET_INCLUDE_DIRECTORIES(_tboxw PRIVATE ${PROJECT_SOURCE_DIR}/fwk/_src
- ${PROJECT_SOURCE_DIR}/tbox/_src
- ${PYTHON_INCLUDE_PATH}
-)
-
-SWIG_LINK_LIBRARIES(tboxw
- tbox fwk ${PYTHON_LIBRARIES}
-)
+TARGET_INCLUDE_DIRECTORIES(tboxw PRIVATE ${PROJECT_SOURCE_DIR}/tbox/_src
+ ${PROJECT_SOURCE_DIR}/fwk/_src
+ ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(tboxw PRIVATE tbox fwk ${PYTHON_LIBRARIES})
INSTALL(FILES "${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/tboxw.py"
DESTINATION ${CMAKE_INSTALL_PREFIX})
-INSTALL(TARGETS _tboxw DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(TARGETS tboxw DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/tbox/_src/tboxw.h b/tbox/_src/tboxw.h
index b9feacdafc52341b820365f648c374ef1b68ced9..3da0640093557e341d636782a8ea0c232b2ca019 100644
--- a/tbox/_src/tboxw.h
+++ b/tbox/_src/tboxw.h
@@ -25,6 +25,7 @@
#include "wCacheQuad4.h"
#include "wCacheTetra4.h"
#include "wCacheTri3.h"
+#include "wDss.h"
#include "wElement.h"
#include "wFct0.h"
#include "wFct1.h"
diff --git a/tbox/_src/tboxw.i b/tbox/_src/tboxw.i
index 36e9d5a8d884fece0dcaf5ea7e41263788d6aaad..135b6a94aeca4184098dd2ad2ef54ffaefd6aa7d 100644
--- a/tbox/_src/tboxw.i
+++ b/tbox/_src/tboxw.i
@@ -97,9 +97,7 @@ namespace std {
%include "wNode.h"
%nodefault tbox::Element; // Element is pure virtual but swig doesn't know it.
-%warnfilter(509); // Shadowed overloaded method MATTYPE. Code compiling and running.
%include "wElement.h"
-%warnfilter(+509);
%include "wTag.h"
%include "wQuad4.h"
%include "wTetra4.h"
@@ -143,7 +141,6 @@ namespace std {
%warnfilter(+509);
%shared_ptr(tbox::MshDeform);
%immutable tbox::MshDeform::msh; // read-only variable (no setter)
-%immutable tbox::MshDeform::nthreads;
%include "wMshDeform.h" // included after std templates because it needs them
%shared_ptr(tbox::MshCrack);
%include "wMshCrack.h"
@@ -156,6 +153,8 @@ namespace std {
%include "wLinearSolver.h"
%shared_ptr(tbox::Pardiso);
%include "wPardiso.h"
+%shared_ptr(tbox::Dss);
+%include "wDss.h"
%shared_ptr(tbox::Mumps);
%include "wMumps.h"
%shared_ptr(tbox::SparseLu);
diff --git a/tbox/src/CMakeLists.txt b/tbox/src/CMakeLists.txt
index c2ccf55085a1e5a9e4382099f7145a9ad2441206..9aedbbc125b0d24212fb4774f66ad4a4923ee6cb 100644
--- a/tbox/src/CMakeLists.txt
+++ b/tbox/src/CMakeLists.txt
@@ -25,12 +25,12 @@ TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${PROJECT_SOURCE_DIR}/tbox/src)
FIND_PACKAGE(TBB REQUIRED)
SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${TBB_CXX_FLAGS_DEBUG}")
TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${TBB_INCLUDE_DIRS})
-IF(TBB_VERSION_INTERFACE LESS 11004)
+IF(TBB_INTERFACE_VERSION LESS 11004)
TARGET_COMPILE_DEFINITIONS(tbox PUBLIC TBB_PREVIEW_GLOBAL_CONTROL)
ENDIF()
TARGET_LINK_LIBRARIES(tbox ${TBB_LIBRARIES})
-# -- MKL/BLAS/LAPACK --
+# -- MKL --
# Do not forget to set your environement variables using: (default path)
# - "C:\Program Files (x86)\Intel\Composer XE\mkl\bin\intel64\mklvars_intel64.bat" (windows, check that VS is in x64 mode)
# - source /opt/intel/mkl/bin/mklvars.sh intel64 (linux)
@@ -45,40 +45,17 @@ IF(USE_MKL)
MESSAGE(STATUS "MKL_LIBRARIES=${MKL_LIBRARIES}")
IF(MKL_INCLUDE_DIRS AND MKL_LIBRARIES)
MESSAGE(STATUS "Found Intel MKL")
- SET(FOUND_MKL TRUE)
- SET(LAPACK_INCLUDE_DIRS ${MKL_INCLUDE_DIRS})
- SET(LAPACK_LIBRARIES ${MKL_LIBRARIES})
+ TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${MKL_INCLUDE_DIRS})
+ TARGET_LINK_LIBRARIES(tbox ${MKL_LIBRARIES})
ELSE()
MESSAGE(FATAL_ERROR "Intel MKL not found!")
ENDIF()
-# Find BLAS
-ELSE()
- IF(${BLA_VENDOR} MATCHES "OpenBlas")
- FIND_LIBRARY(LAPACK_LIBRARIES openblas)
- ELSE()
- FIND_PACKAGE(LAPACK REQUIRED) #FIND_LAPACK calls FIND_BLAS
- ENDIF()
- IF(NOT LAPACK_LIBRARIES)
- MESSAGE(FATAL_ERROR "BLAS/LAPACK not found!")
- ENDIF()
- MESSAGE(STATUS "LAPACK_LIBRARIES=${LAPACK_LIBRARIES}")
- MESSAGE(STATUS "LAPACK_LINKER_FLAGS=${LAPACK_LINKER_FLAGS}")
- MESSAGE(STATUS "BLA_VENDOR=${BLA_VENDOR}")
- MESSAGE(STATUS "BLA_STATIC=${BLA_STATIC}")
ENDIF()
# -- Eigen --
FIND_PACKAGE(EIGEN 3.3.4 REQUIRED)
-TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${EIGEN_INCLUDE_DIRS} ${LAPACK_INCLUDE_DIRS})
+TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${EIGEN_INCLUDE_DIRS})
TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_DONT_PARALLELIZE)
-IF(FOUND_MKL)
- MESSAGE(STATUS "Linking Eigen with MKL")
- TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_USE_MKL_ALL)
-ELSE()
- MESSAGE(STATUS "Linking Eigen with BLAS ${BLA_VENDOR}")
- TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_USE_BLAS)
-ENDIF()
-TARGET_LINK_LIBRARIES(tbox ${LAPACK_LIBRARIES})
# -- MUMPS --
IF (USE_MUMPS)
diff --git a/tbox/src/tbox.h b/tbox/src/tbox.h
index b6814d4c6444cc2afffd02c09a043b438db500b7..d13650219ec96bdffc6662d3759f9a2c19ce4ebd 100644
--- a/tbox/src/tbox.h
+++ b/tbox/src/tbox.h
@@ -99,6 +99,7 @@ class Linesearch;
// Linear solvers
class LinearSolver;
class Pardiso;
+class Dss;
class Mumps;
class SparseLu;
class Gmres;
diff --git a/tbox/src/wDss.cpp b/tbox/src/wDss.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3aaf5235356874006cedd45b5d30a5b8c50df7c5
--- /dev/null
+++ b/tbox/src/wDss.cpp
@@ -0,0 +1,100 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "wDss.h"
+#ifdef USE_MKL
+
+using namespace tbox;
+
+Dss::Dss() : LinearSolver()
+{
+ handle = NULL;
+ MKL_INT opt = MKL_DSS_MSG_LVL_INFO + MKL_DSS_TERM_LVL_FATAL + MKL_DSS_ZERO_BASED_INDEXING; // continue exec until non-recoverable (fatal) error
+ print(dss_create(handle, opt));
+}
+Dss::~Dss()
+{
+ if (handle)
+ {
+ MKL_INT opt = MKL_DSS_DEFAULTS;
+ print(dss_delete(handle, opt));
+ handle = NULL;
+ }
+}
+
+void Dss::analyze(Eigen::SparseMatrix<double> const &A)
+{
+ // Check that matrix A is square
+ if (A.rows() != A.cols())
+ throw std::runtime_error("tbox::Dss: Matrix A is not square!\n");
+ // Convert to CRS
+ Eigen::SparseMatrix<double, Eigen::RowMajor> AA = A;
+ if (!AA.isCompressed())
+ AA.makeCompressed();
+
+ // Build data structure
+ MKL_INT opt = MKL_DSS_NON_SYMMETRIC; // consider MKL_DSS_SYMMETRIC and MKL_DSS_SYMMETRIC_STRUCTURE
+ MKL_INT n = static_cast<int>(AA.rows()); // matrix size (number of rows and columns)
+ MKL_INT nz = static_cast<int>(AA.nonZeros()); // number of non-zero entries
+ const MKL_INT *rStart = AA.outerIndexPtr(); // starting index of non-zero entry for a row
+ const MKL_INT *cIndex = AA.innerIndexPtr(); // index of column of non-zero entry
+ print(dss_define_structure(handle, opt, rStart, n, n, cIndex, nz));
+
+ // Define reordering
+ opt = MKL_DSS_AUTO_ORDER;
+ print(dss_reorder(handle, opt, 0));
+}
+void Dss::factorize(Eigen::SparseMatrix<double> const &A)
+{
+ Eigen::SparseMatrix<double, Eigen::RowMajor> AA = A;
+ MKL_INT opt = MKL_DSS_INDEFINITE;
+ const double *val = AA.valuePtr();
+ print(dss_factor_real(handle, opt, val));
+}
+void Dss::solve(Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x)
+{
+ MKL_INT opt = MKL_DSS_DEFAULTS;
+ MKL_INT n = 1; // number of rhs
+ const double *rhs = b.data();
+ double *sol = x.data();
+ print(dss_solve_real(handle, opt, rhs, n, sol));
+}
+void Dss::compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x)
+{
+ this->analyze(A);
+ this->factorize(A);
+ this->solve(b, x);
+}
+
+/**
+ * @brief Print status code
+ */
+void Dss::print(int cod)
+{
+ if (cod != MKL_DSS_SUCCESS)
+ {
+ std::stringstream err;
+ err << "tbox::Dss: DSS returned code " << cod << "\n";
+ throw std::runtime_error(err.str());
+ }
+}
+
+void Dss::write(std::ostream &out) const
+{
+ out << "Intel MKL DSS (PARDISO)" << std::endl;
+}
+
+#endif // USE_MKL
diff --git a/tbox/src/wDss.h b/tbox/src/wDss.h
new file mode 100644
index 0000000000000000000000000000000000000000..da42e4ac37db46fd320c6bef2e4a8814db7278d2
--- /dev/null
+++ b/tbox/src/wDss.h
@@ -0,0 +1,62 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef WDSS_H
+#define WDSS_H
+
+#include "tbox.h"
+#ifdef USE_MKL
+
+#include "wLinearSolver.h"
+#include <mkl_dss.h>
+#include <Eigen/Sparse>
+#include <limits>
+
+namespace tbox
+{
+/**
+ * @brief Intel DSS (PARDISO) interface
+ * @authors Adrien Crovato
+ */
+class TBOX_API Dss : public LinearSolver
+{
+ void *handle; // DSS handle
+
+ void print(int cod);
+
+public:
+ Dss();
+ ~Dss();
+ virtual LSolType type() const override { return LSolType::DSS; }
+
+#ifndef SWIG
+ virtual void analyze(Eigen::SparseMatrix<double> const &A) override;
+ virtual void factorize(Eigen::SparseMatrix<double> const &A) override;
+ virtual void solve(Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x) override;
+ virtual void compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x) override;
+
+ virtual double getError() override { return std::numeric_limits<double>::epsilon(); }
+ virtual int getIterations() override { return 1; }
+
+ virtual void write(std::ostream &out) const override;
+#endif
+};
+
+} // namespace tbox
+
+#endif // USE_MKL
+
+#endif // WDSS_H
\ No newline at end of file
diff --git a/tbox/src/wElement.cpp b/tbox/src/wElement.cpp
index f2755ff8077b5d5bdde23d0ec6d4f60e3b3dddba..38b1c789a9dd1102ce549ce6e6924bdcdbe1c32e 100644
--- a/tbox/src/wElement.cpp
+++ b/tbox/src/wElement.cpp
@@ -26,26 +26,26 @@ namespace tbox
{
TBOX_API std::ostream &
-operator<<(std::ostream &out, ELTYPE const &obj)
+operator<<(std::ostream &out, ElType const &obj)
{
switch (obj)
{
- case ELTYPE::LINE2:
+ case ElType::LINE2:
out << "LINE2";
break;
- case ELTYPE::TRI3:
+ case ElType::TRI3:
out << "TRI3";
break;
- case ELTYPE::QUAD4:
+ case ElType::QUAD4:
out << "QUAD4";
break;
- case ELTYPE::TETRA4:
+ case ElType::TETRA4:
out << "TETRA4";
break;
- case ELTYPE::HEX8:
+ case ElType::HEX8:
out << "HEX8";
break;
- case ELTYPE::POINT1:
+ case ElType::POINT1:
out << "POINT1";
break;
default:
diff --git a/tbox/src/wElement.h b/tbox/src/wElement.h
index 3c02267a31885b6d3f8bcc961503387628244f38..000c0fc8f7b1e8b4f0fd9a16b1833616b7a4293b 100644
--- a/tbox/src/wElement.h
+++ b/tbox/src/wElement.h
@@ -29,7 +29,7 @@ namespace tbox
* @brief element types (from gmsh)
* @authors Romain Boman
*/
-enum class ELTYPE
+enum class ElType
{
UNKNOWN = 0,
LINE2 = 1,
@@ -40,7 +40,7 @@ enum class ELTYPE
POINT1 = 15
};
#ifndef SWIG
-TBOX_API std::ostream &operator<<(std::ostream &out, ELTYPE const &obj);
+TBOX_API std::ostream &operator<<(std::ostream &out, ElType const &obj);
#endif
/**
@@ -81,9 +81,9 @@ public:
Element(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
virtual ~Element() {}
#endif
- virtual ELTYPE type() const
+ virtual ElType type() const
{
- return ELTYPE::UNKNOWN;
+ return ElType::UNKNOWN;
}
#ifndef SWIG
diff --git a/tbox/src/wGmres.cpp b/tbox/src/wGmres.cpp
index b708c6e54fb045932d808cc5ddcbd9ef39158fa2..40a7db3c275cf368865dd5bc15e948e7f8a3a378 100644
--- a/tbox/src/wGmres.cpp
+++ b/tbox/src/wGmres.cpp
@@ -17,14 +17,14 @@
#include "wGmres.h"
using namespace tbox;
-Gmres::Gmres() : LinearSolver()
+Gmres::Gmres(int _fill, double _tsh, int _rst, double _tol, int _mit) : LinearSolver(), fill(_fill), tsh(_tsh), rst(_rst), tol(_tol), mit(_mit)
{
- // set the preconditioner and solver parameters to "good" starting values
- solver.preconditioner().setFillfactor(1);
- solver.preconditioner().setDroptol(1e-6);
- solver.set_restart(50);
- solver.setTolerance(1e-8);
- solver.setMaxIterations(1000);
+ // set the preconditioner and solver parameters
+ solver.preconditioner().setFillfactor(fill);
+ solver.preconditioner().setDroptol(tsh);
+ solver.set_restart(rst);
+ solver.setTolerance(tol);
+ solver.setMaxIterations(mit);
}
void Gmres::analyze(Eigen::SparseMatrix<double> const &A)
@@ -45,7 +45,13 @@ void Gmres::compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::Vect
x = solver.solve(b);
}
+void Gmres::setTolerance(double t)
+{
+ tol = t;
+ solver.setTolerance(t);
+}
+
void Gmres::write(std::ostream &out) const
{
- out << "GMRES (Eigen)" << std::endl;
+ out << "Eigen GMRES(" << rst << "," << static_cast<int>(log10(tol)) << ")/ILUT(" << fill << "," << static_cast<int>(log10(tsh)) << ")" << std::endl;
}
diff --git a/tbox/src/wGmres.h b/tbox/src/wGmres.h
index b0a4a44ca14328533837ca4d3ace03afc91b33fe..928bf34fdceb9fe4854d675350c158ee9e5015d6 100644
--- a/tbox/src/wGmres.h
+++ b/tbox/src/wGmres.h
@@ -30,11 +30,18 @@ namespace tbox
*/
class TBOX_API Gmres : public LinearSolver
{
+ int fill; // fill-in factor for ILUT
+ double tsh; // drop tolerance threshold for ILUT
+ int rst; // number of iterations before restart for GMRES
+ double tol; // relative tolerance for GMRES
+ int mit; // maximum number of iterations for GMRES
Eigen::GMRES<Eigen::SparseMatrix<double>, Eigen::IncompleteLUT<double>> solver;
public:
- Gmres();
+ Gmres(int _fill = 1, double _tsh = 1e-6, int _rst = 30, double _tol = 1e-8, int _mit = 1000);
+ Gmres(const Gmres &gmres) : Gmres(gmres.fill, gmres.tsh, gmres.rst, gmres.tol, gmres.mit) {}
virtual ~Gmres() {}
+ virtual LSolType type() const override { return LSolType::GMRES_ILUT; }
#ifndef SWIG
virtual void analyze(Eigen::SparseMatrix<double> const &A) override;
@@ -43,19 +50,12 @@ public:
virtual void compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b, Eigen::Map<Eigen::VectorXd> &x) override;
virtual double getError() override { return solver.error(); }
- virtual int getIterations() override { return solver.iterations(); }
+ virtual int getIterations() override { return static_cast<int>(solver.iterations()); }
+
+ void setTolerance(double t);
virtual void write(std::ostream &out) const override;
#endif
-
- void setFillFactor(int f)
- {
- solver.preconditioner().setFillfactor(f);
- }
- void setDropTol(double t) { solver.preconditioner().setDroptol(t); }
- void setTolerance(double t) { solver.setTolerance(t); }
- void setRestart(int r) { solver.set_restart(r); }
- void setMaxIterations(int n) { solver.setMaxIterations(n); }
};
} // namespace tbox
diff --git a/tbox/src/wGmshExport.h b/tbox/src/wGmshExport.h
index eefcf4470d05730ceeaa7aba8552ae93c06c15bd..dc74bf3ebb68693991f9e7f82bba716344b45f9c 100644
--- a/tbox/src/wGmshExport.h
+++ b/tbox/src/wGmshExport.h
@@ -34,6 +34,7 @@ class TBOX_API GmshExport : public MshExport
public:
GmshExport(std::shared_ptr<MshData> _msh);
virtual ~GmshExport();
+ virtual WrtType type() const override { return WrtType::GMSH; }
virtual void save(std::string const &fname) const override;
#ifndef SWIG
@@ -45,4 +46,4 @@ public:
} // namespace tbox
-#endif //WGMSHUTILS_H
+#endif // WGMSHEXPORT_H
diff --git a/tbox/src/wGmshImport.cpp b/tbox/src/wGmshImport.cpp
index 3ea091dcafff117000b8284cbb216c7821ff2c80..4a03444fa641e20bca9cbc0cdfdc8308846b875b 100644
--- a/tbox/src/wGmshImport.cpp
+++ b/tbox/src/wGmshImport.cpp
@@ -256,7 +256,7 @@ void GmshImport::readMshElements(FILE *file, int myrank)
belongsToPart = true;
}
}
- if (static_cast<ELTYPE>(type) == ELTYPE::POINT1)
+ if (static_cast<ElType>(type) == ElType::POINT1)
belongsToPart = true;
bool skipline = true;
@@ -265,32 +265,32 @@ void GmshImport::readMshElements(FILE *file, int myrank)
{
skipline = false;
std::vector<Node *> nods;
- if (static_cast<ELTYPE>(type) == ELTYPE::HEX8) // 8-node hexa
+ if (static_cast<ElType>(type) == ElType::HEX8) // 8-node hexa
{
readElNods(file, 8, nods);
msh->elems.push_back(new Hex8(no, ptag, etag, prts, nods));
}
- else if (static_cast<ELTYPE>(type) == ELTYPE::QUAD4) // 4-node quad
+ else if (static_cast<ElType>(type) == ElType::QUAD4) // 4-node quad
{
readElNods(file, 4, nods);
msh->elems.push_back(new Quad4(no, ptag, etag, prts, nods));
}
- else if (static_cast<ELTYPE>(type) == ELTYPE::TETRA4) // 4-node tetra
+ else if (static_cast<ElType>(type) == ElType::TETRA4) // 4-node tetra
{
readElNods(file, 4, nods);
msh->elems.push_back(new Tetra4(no, ptag, etag, prts, nods));
}
- else if (static_cast<ELTYPE>(type) == ELTYPE::LINE2) // 2-node line
+ else if (static_cast<ElType>(type) == ElType::LINE2) // 2-node line
{
readElNods(file, 2, nods);
msh->elems.push_back(new Line2(no, ptag, etag, prts, nods));
}
- else if (static_cast<ELTYPE>(type) == ELTYPE::POINT1) // 1-node point
+ else if (static_cast<ElType>(type) == ElType::POINT1) // 1-node point
{
readElNods(file, 1, nods);
msh->elems.push_back(new Point1(no, ptag, etag, prts, nods));
}
- else if (static_cast<ELTYPE>(type) == ELTYPE::TRI3) // 3-node triangle
+ else if (static_cast<ElType>(type) == ElType::TRI3) // 3-node triangle
{
readElNods(file, 3, nods);
msh->elems.push_back(new Tri3(no, ptag, etag, prts, nods));
diff --git a/tbox/src/wGroup.cpp b/tbox/src/wGroup.cpp
index c26634fd859ab247d3e102ea628a242f8ff84203..b271543f2441c8da0d23f96aea35a90753eaa434 100644
--- a/tbox/src/wGroup.cpp
+++ b/tbox/src/wGroup.cpp
@@ -28,7 +28,7 @@ Group::Group(std::shared_ptr<MshData> _msh, int no) : wSharedObject(), msh(_msh)
if (it == msh->ptags.end())
{
std::stringstream out;
- out << "Physical Group #" << no << " not found among physical tags";
+ out << "Group: Physical Group #" << no << " not found among physical tags";
throw std::out_of_range(out.str());
}
tag = it->second;
@@ -40,8 +40,8 @@ Group::Group(std::shared_ptr<MshData> _msh, std::string const &name) : wSharedOb
if (it == msh->ntags.end())
{
std::stringstream out;
- out << "Physical Group \"" << name << "\" not found among physical tags\n";
- out << "[info] available tags:\n\t";
+ out << "Group: Physical Group \"" << name << "\" not found among physical tags\n";
+ out << "[info] available tags: ";
for (auto it = msh->ntags.begin(); it != msh->ntags.end(); ++it)
{
auto itn(it);
diff --git a/tbox/src/wHex8.h b/tbox/src/wHex8.h
index 391b73dce4b71ae7b5866f7d2b26dfc8bf1fe19b..fb3e427b5cc27efdd7e51f709c7cbffe75b9083e 100644
--- a/tbox/src/wHex8.h
+++ b/tbox/src/wHex8.h
@@ -43,7 +43,7 @@ protected:
public:
Hex8(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
virtual ~Hex8() {}
- virtual ELTYPE type() const override { return ELTYPE::HEX8; }
+ virtual ElType type() const override { return ElType::HEX8; }
#ifndef SWIG
virtual void initValues(bool isvol, int ordr = 2) override;
diff --git a/tbox/src/wLine2.h b/tbox/src/wLine2.h
index 6b076608d104fc2ffe5877ee05ed9a69c918d68a..d87928f5898f76558d3161d4493790a26d7343b6 100644
--- a/tbox/src/wLine2.h
+++ b/tbox/src/wLine2.h
@@ -49,7 +49,7 @@ protected:
public:
Line2(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
virtual ~Line2() {}
- virtual ELTYPE type() const override { return ELTYPE::LINE2; }
+ virtual ElType type() const override { return ElType::LINE2; }
#ifndef SWIG
virtual void initValues(bool isvol, int ordr = 2) override;
diff --git a/tbox/src/wLinearSolver.h b/tbox/src/wLinearSolver.h
index 3e93813a4db7b1899b03c5fd6cfb80b81c44bfd6..b717b8c8c538a943dd73439da5df0e56aeda50c3 100644
--- a/tbox/src/wLinearSolver.h
+++ b/tbox/src/wLinearSolver.h
@@ -24,6 +24,20 @@
namespace tbox
{
+/**
+ * @brief Solver type
+ * @authors Adrien Crovato
+ */
+enum class LSolType
+{
+ UNDEFINED = 0,
+ SPARSE_LU = 1,
+ GMRES_ILUT = 2,
+ PARDISO = 3,
+ DSS = 4,
+ MUMPS = 5
+};
+
/**
* @brief Base interface class for linear solvers
* @authors Adrien Crovato
@@ -33,6 +47,7 @@ class TBOX_API LinearSolver : public fwk::wSharedObject
public:
LinearSolver() {}
virtual ~LinearSolver() {}
+ virtual LSolType type() const { return LSolType::UNDEFINED; }
#ifndef SWIG
virtual void analyze(Eigen::SparseMatrix<double> const &A);
diff --git a/tbox/src/wLinesearch.cpp b/tbox/src/wLinesearch.cpp
index 3bf8e33032bce9b003ac5c70f189202e267541b5..8254d77f9e6780993fe0a79158efcf55029c1d70 100644
--- a/tbox/src/wLinesearch.cpp
+++ b/tbox/src/wLinesearch.cpp
@@ -111,7 +111,7 @@ double FleuryLS::run()
double phit = fct.eval(h);
fevalIt = 2;
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Starting Backeting \n";
std::cout << "\t" << std::setprecision(16) << "phi1(" << 0 << ") = " << phi1 << '\n';
@@ -128,14 +128,14 @@ double FleuryLS::run()
{
if (h > hU)
{
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Positive slope detected (in ascending bracketing #1)!\n";
}
if (revAllowed)
{
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Reverse Direction of Line Search (in ascending bracketing #1)!"
<< "\n";
@@ -146,7 +146,7 @@ double FleuryLS::run()
}
else
{
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Error in ascending bracketing #1:\n"
<< std::setprecision(16) << "a1=" << 0 << " a2=" << h << " a3 =" << 2.0 * h << '\n'
@@ -164,7 +164,7 @@ double FleuryLS::run()
phit = fct.eval(2 * h);
fevalIt++;
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Ascending bracketing iteration # " << nB << ": \n";
std::cout << "\t"
@@ -185,14 +185,14 @@ double FleuryLS::run()
{
if (h < hL)
{
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Positive slope detected (in descending bracketing #2)!\n";
}
if (revAllowed)
{
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Reverse Direction of Line Search (in descending bracketing #2) !"
<< "\n";
@@ -203,7 +203,7 @@ double FleuryLS::run()
}
else
{
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Error in descending bracketing #2:\n"
<< std::setprecision(16) << "a1=" << 0 << " a2=" << h << " a3 =" << 2.0 * h << '\n'
@@ -220,7 +220,7 @@ double FleuryLS::run()
phi3 = phit;
phit = fct.eval(h / 2);
fevalIt++;
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Descending bracketing iteration # " << nB << ": \n";
std::cout << "\t"
@@ -241,7 +241,7 @@ double FleuryLS::run()
double a3 = 2 * h;
double a4 = h * (4 * phi2 - 3 * phi1 - phi3) / (4 * phi2 - 2 * phi1 - 2 * phi3);
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Backeting succeeds in " << nB << " iterations : \n";
std::cout << "\t" << std::setprecision(16) << "phi1(" << a1 << ") = " << phi1 << '\n';
@@ -276,7 +276,7 @@ double FleuryLS::run()
throw std::runtime_error(msg.str());
}
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Line search iteration # " << nLS << ": \n";
std::cout << "\t" << std::setprecision(16) << "phi1(" << a1 << ") = " << phi1 << '\n';
@@ -319,7 +319,7 @@ double FleuryLS::run()
A = (a2 - a3) * phi1 + (a3 - a1) * phi2 + (a1 - a2) * phi3;
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "\t" << std::setprecision(16) << "|A| = " << fabs(A) << "><" << tol << '\n';
}
@@ -337,7 +337,7 @@ double FleuryLS::run()
(A);
}
- if (verbose > 2)
+ if (verbose > 3)
{
std::cout << "Line search succeeds in " << nLS << " iterations : \n";
std::cout << "\t" << std::setprecision(16) << "|A| = " << fabs(A) << '\n';
diff --git a/tbox/src/wMshConvert.cpp b/tbox/src/wMshConvert.cpp
index 013fff0c07c8539820a1d14ae2fd54138aafa852..dd48f05554736b46383083f81a0a29fa934adb95 100644
--- a/tbox/src/wMshConvert.cpp
+++ b/tbox/src/wMshConvert.cpp
@@ -22,6 +22,7 @@
#include <fstream>
#include <set>
#include <algorithm>
+#include <limits>
using namespace tbox;
diff --git a/tbox/src/wMshCrack.cpp b/tbox/src/wMshCrack.cpp
index 9e11b8a8f62a14dbb8d5b511bc299795a40d08bd..6ab84c00b74ff50cec732b8b0cf502bb039f68a0 100644
--- a/tbox/src/wMshCrack.cpp
+++ b/tbox/src/wMshCrack.cpp
@@ -28,7 +28,7 @@ using namespace tbox;
MshCrack::MshCrack(std::shared_ptr<MshData> _msh, int _nDim) : wSharedObject(), msh(_msh), nDim(_nDim)
{
- hasCrack = true; // so that we prevent automatic data structure creation if addCrack() is not called
+ hasCrack = true; // so that we prevent automatic data structure creation if setCrack() is not called
crckGrp = nullptr;
exclGrp = nullptr;
}
@@ -37,58 +37,84 @@ MshCrack::~MshCrack()
// Delete temporary groups
delete crckGrp;
delete exclGrp;
- for (size_t i = 0; i < bndGrps.size(); ++i)
- delete bndGrps[i];
- for (size_t i = 0; i < bndGrp.size(); ++i)
- delete bndGrp[i];
+ for (size_t i = 0; i < sBndGrps.size(); ++i)
+ delete sBndGrps[i];
+ for (size_t i = 0; i < sBndGrp.size(); ++i)
+ delete sBndGrp[i];
+ for (size_t i = 0; i < vBndGrps.size(); ++i)
+ delete vBndGrps[i];
+ for (size_t i = 0; i < vBndGrp.size(); ++i)
+ delete vBndGrp[i];
std::cout << "~MshCrack()\n";
}
/**
* @brief Create temporary group for crack
*/
-void MshCrack::setCrack(std::string const &crckName)
+void MshCrack::setCrack(std::string const &name)
{
// Check if a crack needs to be made
- auto it = msh->ntags.find(crckName + '_');
+ auto it = msh->ntags.find(name + '_');
if (it != msh->ntags.end())
- std::cout << "Groups " << crckName << " and " << crckName << "_ are already existing. MshCrack will not be run!\n";
+ std::cout << "Groups " << name << " and " << name << "_ are already existing. MshCrack will not be run\n";
else
{
hasCrack = false;
- crckGrp = new Group(msh, crckName);
+ crckGrp = new Group(msh, name);
+ if (crckGrp->tag->dim != nDim - 1)
+ {
+ std::stringstream err;
+ err << "MshCrack::setCrack: cannot create a crack of dimension " << crckGrp->tag->dim << ", dimension must be " << nDim - 1 << "!\n";
+ throw std::runtime_error(err.str());
+ }
}
}
/**
* @brief Create temporary group for excluded nodes on crack
*/
-void MshCrack::setExcluded(std::string const &exclName)
+void MshCrack::setExcluded(std::string const &name)
{
if (!hasCrack)
- exclGrp = new Group(msh, exclName);
+ exclGrp = new Group(msh, name);
}
/**
* @brief Create temporary groups for crack boundaries
*/
-void MshCrack::addBoundaries(std::vector<std::string> const &bndName)
+void MshCrack::addBoundary(std::string const &name)
{
if (!hasCrack)
{
- for (auto name : bndName)
+ // Check if the boundaries on both sides of the crack are explicitely provided in two different tags
+ // if so, store them
+ try
{
- // Check if the boundaries are explicitely separated...
- try
+ Groups *gs = new Groups(msh, std::vector<std::string>{name, name + '_'});
+ if (gs->groups[0]->tag->dim == nDim)
+ vBndGrps.push_back(gs);
+ else if (gs->groups[0]->tag->dim == nDim - 1)
+ sBndGrps.push_back(gs);
+ else
{
- Groups *gs = new Groups(msh, std::vector<std::string>{name, name + '_'});
- bndGrps.push_back(gs);
+ std::stringstream err;
+ err << "MshCrack::addBoundary: cannot add a boundary of dimension " << gs->groups[0]->tag->dim << ", dimension must be either " << nDim << " or " << nDim - 1 << "!\n";
+ throw std::runtime_error(err.str());
}
- // ... if not, an algorithm will later try to to separate the elements according to the side of the crack they are touching
- catch (const std::out_of_range &)
+ }
+ // else, the boundaries are contained in a single tag, an algorithm will later try to separate the elements according to the side of the crack they are on
+ catch (const std::out_of_range &)
+ {
+ Group *g = new Group(msh, name);
+ if (g->tag->dim == nDim)
+ vBndGrp.push_back(g);
+ else if (g->tag->dim == nDim - 1)
+ sBndGrp.push_back(g);
+ else
{
- Group *g = new Group(msh, name);
- bndGrp.push_back(g);
+ std::stringstream err;
+ err << "MshCrack::addBoundary: cannot add a boundary of dimension " << g->tag->dim << ", dimension must be either " << nDim << " or " << nDim - 1 << "!\n";
+ throw std::runtime_error(err.str());
}
}
}
@@ -100,15 +126,14 @@ void MshCrack::addBoundaries(std::vector<std::string> const &bndName)
void MshCrack::run()
{
if (hasCrack)
- std::cout << "MshCrack not run!\n";
+ std::cout << "MshCrack not run\n";
else
{
// Check that crack is bounded
- if (bndGrps.empty() && bndGrp.empty())
- throw std::runtime_error("MshCrack::run crack has no boundary (use MshCrack::addBoundary()!\n");
+ if ((sBndGrps.empty() && sBndGrp.empty()) || (vBndGrps.empty() && vBndGrp.empty()))
+ throw std::runtime_error("MshCrack::run, crack has no (surface and/or volume) boundary, use MshCrack::addBoundary()!\n");
// Open the crack and modify data structure
openCrack();
- modifyBoundaries_();
modifyBoundaries();
std::cout << "Crack added on " << *(crckGrp->tag) << "(" << crckGrp->tag->name << "_ created)" << std::endl;
}
@@ -177,14 +202,14 @@ void MshCrack::openCrack()
nods[j] = e->nodes[e->nodes.size() - j - 1];
}
}
- if (e->type() == ELTYPE::LINE2 && nDim == 2)
+ if (e->type() == ElType::LINE2 && nDim == 2)
msh->elems.push_back(new Line2(static_cast<int>(i), tagp, tage, e->parts, nods));
- else if (e->type() == ELTYPE::TRI3 && nDim == 3)
+ else if (e->type() == ElType::TRI3 && nDim == 3)
msh->elems.push_back(new Tri3(static_cast<int>(i), tagp, tage, e->parts, nods));
else
{
std::stringstream err;
- err << "Msh::Crack: could not create element of type " << e->type() << " for dimension " << nDim << "!\n";
+ err << "MshCrack::run: could not create element of type " << e->type() << " for dimension " << nDim << "!\n";
throw std::runtime_error(err.str());
}
i++;
@@ -192,15 +217,123 @@ void MshCrack::openCrack()
}
/**
- * @brief Swap the elements nodes on one side of the crack and merge physical tags
+ * @brief Modify the elements on the negative side of the crack
+ *
+ * 1) Identify boundary surface elements on negative side of crack, and store their nodes
+ * - easy case: elements are contained in different tags
+ * - hard case: elements are contained in a single tag, they are identified by comparing the vector joining them to the crack to the normal vector of the crack (NOT ROBUST)
+ * 2) Identify boundary volume elements on negative side of crack
+ * - easy case: elements are contained in different tags
+ * - hard case: elements are contained in a single tag, they are identified by
+ * - either checking if one of their nodes belong to the surface
+ * - or, comparing the vector joining them to the crack to the normal vector of the crack
+ * 3) Swap the nodes of (surface and volume) elements located on the negative side of the crack
+ * 4) Merge the tags identifying the (surface and volume) elements on the negative side into the tags identifying the positive side
*/
-void MshCrack::modifyBoundaries_()
+void MshCrack::modifyBoundaries()
{
- for (auto gs : bndGrps)
+ // Store crack nodes, and compute CGs and normals if sides are contained in a single tag
+ std::unordered_set<Node *> crkNods, crkNodsAll;
+ std::vector<Eigen::Vector3d> crkCG(crckGrp->tag->elems.size()), crkN(crckGrp->tag->elems.size());
+ if (!sBndGrp.empty() || !vBndGrp.empty())
+ {
+ // nodes of crack except excluded
+ for (auto p : nodMap)
+ crkNods.insert(p.first);
+ // all nodes of crack
+ for (auto e : crckGrp->tag->elems)
+ crkNodsAll.insert(e->nodes.begin(), e->nodes.end());
+ // cg and unit normal vector
+ for (size_t i = 0; i < crckGrp->tag->elems.size(); ++i)
+ {
+ crckGrp->tag->elems[i]->computeCg();
+ crkCG[i] = crckGrp->tag->elems[i]->cg;
+ crckGrp->tag->elems[i]->computeNormal();
+ crkN[i] = crckGrp->tag->elems[i]->normal;
+ }
+ }
+
+ // Identify surface elements on negative side, and store their nodes (to find volume elements if they are contained in a single tag)
+ std::vector<Element *> _surEl;
+ std::unordered_set<Node *> _surNd;
+ // Each side is contained in a different tag
+ for (auto gs : sBndGrps)
{
- // Modify element nodes touching the crack
for (auto e : gs->groups[1]->tag->elems)
- swapNodes(e);
+ {
+ _surEl.push_back(e);
+ // only keep nodes that are not on crack (since they are on positive side)
+ for (auto n : e->nodes)
+ if (!crkNodsAll.count(n))
+ _surNd.insert(n);
+ }
+ }
+ // Both sides are contained in a single tag
+ for (auto g : sBndGrp)
+ {
+ // find elements having at least one node on the crack, and store if on negative side
+ std::unordered_set<Element *> onCrk = findTouching(g->tag->elems, crkNods);
+ for (auto e : onCrk)
+ {
+ // find mininal distance bewteen CGs
+ e->computeCg();
+ size_t idx = findClosest(crkCG, e->cg);
+ // if dot product between crack-to-volume vector and crack normal is negative, element is on negative side
+ if ((e->cg - crkCG[idx]).dot(crkN[idx]) < 0)
+ {
+ _surEl.push_back(e);
+ // only keep nodes that are not on crack (since they are on positive side)
+ for (auto n : e->nodes)
+ if (!crkNodsAll.count(n))
+ _surNd.insert(n);
+ }
+ }
+ }
+
+ // Identify volume elements touching surface elements on negative side
+ std::vector<Element *> _volEl;
+ // Each side is contained in a different tag
+ for (auto gs : vBndGrps)
+ for (auto e : gs->groups[1]->tag->elems)
+ _volEl.push_back(e);
+ // Both sides are contained in a single tag
+ for (auto g : vBndGrp)
+ {
+ // find elements having at least one node on the crack
+ std::unordered_set<Element *> onCrkOnly = findTouching(g->tag->elems, crkNods);
+ // find elements having at least one node on the negative side of the surface
+ std::unordered_set<Element *> onSurf = findTouching(std::vector<Element *>(onCrkOnly.begin(), onCrkOnly.end()), _surNd);
+ // remove these elements from the first set, as they are necessarily on the negative side
+ _volEl.insert(_volEl.end(), onSurf.begin(), onSurf.end());
+ for (auto e : onSurf)
+ onCrkOnly.erase(e);
+ // check if remaining elements are on the negative side
+ for (auto e : onCrkOnly)
+ {
+ // find mininal distance bewteen CGs
+ e->computeCg();
+ size_t idx = findClosest(crkCG, e->cg);
+ // if dot product between crack-to-volume vector and crack normal is negative, element is on negative side
+ if ((e->cg - crkCG[idx]).dot(crkN[idx]) < 0)
+ _volEl.push_back(e);
+ }
+ }
+
+ // Swap nodes of surface and volume elements on negative side
+ swapNodes(_surEl);
+ swapNodes(_volEl);
+ // Merge tag identifying negative side into positive side tag
+ mergeTags(sBndGrps);
+ mergeTags(vBndGrps);
+}
+
+/**
+ * @brief Merge physical tags and modify data structure accordingly
+ */
+void MshCrack::mergeTags(std::vector<Groups *> &grps)
+{
+ for (auto gs : grps)
+ {
// Merge tags
Tag *refTag = gs->groups[0]->tag;
Tag *oldTag = gs->groups[1]->tag;
@@ -217,74 +350,58 @@ void MshCrack::modifyBoundaries_()
}
/**
- * @brief Detect the elements on one side of the crack and swap their nodes
+ * @brief Find elements touching a surface (ie, that have at least a node in the provided set)
*/
-void MshCrack::modifyBoundaries()
+std::unordered_set<Element *> MshCrack::findTouching(std::vector<Element *> const &elms, std::unordered_set<Node *> const &nods)
{
- // compute crack CGs and normals
- std::vector<Eigen::Vector3d> crkCG(crckGrp->tag->elems.size()), crkN(crckGrp->tag->elems.size());
- for (size_t i = 0; i < crckGrp->tag->elems.size(); ++i)
+ std::unordered_set<Element *> found;
+ for (auto e : elms)
{
- crckGrp->tag->elems[i]->computeCg();
- crkCG[i] = crckGrp->tag->elems[i]->cg;
- crckGrp->tag->elems[i]->computeNormal();
- crkN[i] = crckGrp->tag->elems[i]->normal;
- }
- // find elements on the opposite side of the crack and swap their nodes
- for (auto g : bndGrp)
- {
- // find elements that have at least one node on the crack and compute their CG
- std::map<Element *, Eigen::Vector3d> onCrk;
- for (auto e : g->tag->elems)
- {
- for (auto n : e->nodes)
- {
- if (nodMap.find(n) != nodMap.end())
- {
- //Eigen::Vector3d cg(0, 0, 0);
- //for (auto n_ : e->nodes)
- // cg += n_->pos;
- //cg /= e->nodes.size();
- //onCrk[e] = cg;
- e->computeCg();
- onCrk[e] = e->cg;
- break;
- }
- }
- }
- // find closest element, check side and swap nodes
- for (auto p : onCrk)
+ for (auto n : e->nodes)
{
- // find mininal distance bewteen CGs
- size_t idx = 0;
- std::vector<Eigen::Vector3d> deltaCG(crkCG.size());
- for (size_t i = 0; i < crkCG.size(); ++i)
+ if (nods.count(n))
{
- deltaCG[i] = p.second - crkCG[i];
- if (deltaCG[i].norm() < deltaCG[idx].norm())
- idx = i;
+ found.insert(e);
+ break;
}
- // check sign of dot product
- if (deltaCG[idx].dot(crkN[idx]) < 0)
- swapNodes(p.first);
}
}
+ return found;
}
/**
- * @brief Replace nodes of elements on opposite side of crack by opposite nodes
+ * @brief Find closest element and returns its index
*/
-void MshCrack::swapNodes(Element *e)
+size_t MshCrack::findClosest(std::vector<Eigen::Vector3d> const &cga, Eigen::Vector3d const &cgb)
{
- for (size_t i = 0; i < e->nodes.size(); ++i)
+ size_t idx = 0;
+ std::vector<Eigen::Vector3d> delta(cga.size());
+ for (size_t i = 0; i < cga.size(); ++i)
{
- try
- {
- e->nodes[i] = nodMap.at(e->nodes[i]);
- }
- catch (const std::out_of_range &)
+ delta[i] = cgb - cga[i];
+ if (delta[i].norm() < delta[idx].norm())
+ idx = i;
+ }
+ return idx;
+}
+
+/**
+ * @brief Replace nodes of boundary elements on negative side of crack by negative crack nodes
+ */
+void MshCrack::swapNodes(std::vector<Element *> &elms)
+{
+ for (auto e : elms)
+ {
+ for (size_t i = 0; i < e->nodes.size(); ++i)
{
- //std::cout << e->nodes[i]->no << "not found in map!\n";
+ try
+ {
+ e->nodes[i] = nodMap.at(e->nodes[i]);
+ }
+ catch (const std::out_of_range &)
+ {
+ // std::cout << e->nodes[i]->no << "not found in map!\n";
+ }
}
}
}
@@ -292,4 +409,4 @@ void MshCrack::swapNodes(Element *e)
void MshCrack::write(std::ostream &out) const
{
out << "MshCrack on " << msh << std::endl;
-}
\ No newline at end of file
+}
diff --git a/tbox/src/wMshCrack.h b/tbox/src/wMshCrack.h
index dcce73d0dd8de221f3ca0fa3c3c22a34ccbb5a48..92338f7c85b1f56aa76763d35794b897b0b3894e 100644
--- a/tbox/src/wMshCrack.h
+++ b/tbox/src/wMshCrack.h
@@ -19,9 +19,10 @@
#include "tbox.h"
#include "wObject.h"
-#include "wMshExport.h"
+#include <Eigen/Dense>
#include <vector>
#include <map>
+#include <unordered_set>
#include <memory>
namespace tbox
@@ -30,9 +31,10 @@ namespace tbox
/**
* @brief Duplicate nodes and elements in a physical group and make mesh consistent
*
- * Similar to plugin "crack" from gmsh but allows more control. Currently coded
+ * Similar to plugin "crack" from gmsh but allows more control. Currently implemented
* only for Line2 and Tri3.
* @authors Adrien Crovato
+ * @todo Automatic identification of boundary surface elements on negative side of crack is not robust (might fail for concave boundaries)
*/
class TBOX_API MshCrack : public fwk::wSharedObject
{
@@ -43,21 +45,25 @@ private:
std::map<Node *, Node *> nodMap; ///< map between nodes on both side of crack
Group *crckGrp; ///< physical group of crack
Group *exclGrp; ///< physical group of nodes not to be duplicated
- std::vector<Groups *> bndGrps; ///< list of pair of physical groups touching crack
- std::vector<Group *> bndGrp; ///< list of physical groups touching crack
+ std::vector<Groups *> sBndGrps; ///< list of pair of surface physical groups touching crack
+ std::vector<Groups *> vBndGrps; ///< list of pair of volume physical groups touching crack
+ std::vector<Group *> sBndGrp; ///< list of surface physical groups touching crack
+ std::vector<Group *> vBndGrp; ///< list of volume physical groups touching crack
void openCrack();
- void modifyBoundaries_();
void modifyBoundaries();
- void swapNodes(Element *e);
+ std::unordered_set<Element *> findTouching(std::vector<Element *> const &elms, std::unordered_set<Node *> const &nods);
+ size_t findClosest(std::vector<Eigen::Vector3d> const &cga, Eigen::Vector3d const &cgb);
+ void mergeTags(std::vector<Groups *> &grps);
+ void swapNodes(std::vector<Element *> &elms);
public:
MshCrack(std::shared_ptr<MshData> _msh, int _nDim);
virtual ~MshCrack();
- void setCrack(std::string const &crckName);
- void addBoundaries(std::vector<std::string> const &bndName);
- void setExcluded(std::string const &exclName);
+ void setCrack(std::string const &name);
+ void setExcluded(std::string const &name);
+ void addBoundary(std::string const &name);
void run();
#ifndef SWIG
@@ -67,4 +73,4 @@ public:
} // namespace tbox
-#endif //WMSHCRACK_H
\ No newline at end of file
+#endif // WMSHCRACK_H
\ No newline at end of file
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index bb740c91b662f7693d089fa5fae9ddb349cb1386..075a4a50e3a6d9438657c5fbaa8e8b5b7a60621e 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -16,6 +16,7 @@
#include "wMshDeform.h"
#include "wMshData.h"
+#include "wLinearSolver.h"
#include "wGroup.h"
#include "wGroups.h"
#include "wTag.h"
@@ -23,19 +24,17 @@
#include "wElement.h"
#include "wCache.h"
#include "wGauss.h"
-#include "wGmres.h"
#include <tbb/global_control.h>
-#include <tbb/parallel_for_each.h>
-#include <tbb/spin_mutex.h>
#include <deque>
using namespace tbox;
MshDeform::MshDeform(std::shared_ptr<MshData> _msh,
- int _nDim, int nthrds) : wSharedObject(), nDim(_nDim),
- field(false), fixed(false), moving(false),
- msh(_msh), nthreads(nthrds)
+ std::shared_ptr<tbox::LinearSolver> _linsol,
+ int _nDim, int nthrds, int vrb) : wSharedObject(),
+ nDim(_nDim), nthreads(nthrds), verbose(vrb),
+ msh(_msh), linsol(_linsol)
{
// Check problem dimension
if (nDim != 2 && nDim != 3)
@@ -49,6 +48,7 @@ MshDeform::MshDeform(std::shared_ptr<MshData> _msh,
// Display configuration
std::cout << "--- Mesh morpher (linear elasticity) ---\n"
+ << "Inner solver: " << *linsol
<< "Number of threads: " << nthreads << "\n"
<< std::endl;
}
@@ -58,37 +58,25 @@ MshDeform::MshDeform(std::shared_ptr<MshData> _msh,
*/
void MshDeform::setField(std::string const &fldName)
{
- // Create temporary field group
- fldGrp = new Group(msh, fldName);
- field = true;
+ fldGrp = new Group(msh, fldName); // temporary field group
}
/**
* @brief Add fixed boundaries
*/
-void MshDeform::addFixed(std::vector<std::string> const &fxdBndName)
+void MshDeform::addFixed(std::string const &fxdBndName)
{
- // Create temporary boundary groups
- for (size_t i = 0; i < fxdBndName.size(); ++i)
- {
- Group *g = new Group(msh, fxdBndName[i]);
- fxdBndGrp.push_back(g);
- }
- fixed = true;
+ Group *g = new Group(msh, fxdBndName); // temporary boundary group
+ fxdBndGrp.push_back(g);
}
/**
* @brief Add moving boundaries
*/
-void MshDeform::addMoving(std::vector<std::string> const &movBndName)
+void MshDeform::addMoving(std::string const &movBndName)
{
- // Create temporary boundary groups
- for (size_t i = 0; i < movBndName.size(); ++i)
- {
- Group *g = new Group(msh, movBndName[i]);
- movBndGrp.push_back(g);
- }
- moving = true;
+ Group *g = new Group(msh, movBndName); // temporary boundary group
+ movBndGrp.push_back(g);
}
/**
@@ -96,23 +84,19 @@ void MshDeform::addMoving(std::vector<std::string> const &movBndName)
*/
void MshDeform::setSymmetry(std::string const &symBndName, size_t _blck)
{
- // Get blocked direction
+ // Get blocked direction and create group
blck = _blck;
if (blck != 0 && blck != 1 && blck != 2)
- throw std::runtime_error("tbox::MshDeform::setSymmetry: blck should be 0, 1 or 2 (for x, y and z direction resp.)!\n");
- // Create temporary boundary groups
- symBndGrp = new Group(msh, symBndName);
+ throw std::runtime_error("MshDeform::setSymmetry: blck should be 0, 1 or 2 (for x, y and z direction resp.)!\n");
+ symBndGrp = new Group(msh, symBndName); // temporary field group
}
/**
* @brief Add internal boundaries
*/
-void MshDeform::addInternal(std::vector<std::string> const &intBndName)
+void MshDeform::addInternal(std::string const &intBndName, std::string const &intBndName_)
{
- // Create temporary internal link
- if (intBndName.size() != 2)
- throw std::runtime_error("MshDeform: input intBndName must be a pair of internal boundaries!\n");
- Groups *gs = new Groups(msh, intBndName);
+ Groups *gs = new Groups(msh, std::vector<std::string>{intBndName, intBndName_}); // temporary internal link group
intBndGrps.push_back(gs);
}
@@ -122,12 +106,12 @@ void MshDeform::addInternal(std::vector<std::string> const &intBndName)
void MshDeform::initialize()
{
// Sanity checks
- if (!field)
- throw std::runtime_error("MshDeform: field has not been found! Use setField() to add a field.\n");
- else if (!fixed)
- throw std::runtime_error("MshDeform: fixed boundaries have not been found! Use addFixed() to add at least one boundary.\n");
- else if (!moving)
- throw std::runtime_error("MshDeform: moving boundaries have not been found! Use addMoving() to add at least one boundary.\n");
+ if (fldGrp == nullptr)
+ throw std::runtime_error("MshDeform::initialize: field has not been found! Use setField() to add a field.\n");
+ if (fxdBndGrp.empty())
+ throw std::runtime_error("MshDeform::initialize: fixed boundaries have not been found! Use addFixed() to add at least one boundary.\n");
+ if (movBndGrp.empty())
+ throw std::runtime_error("MshDeform::initialize: moving boundaries have not been found! Use addMoving() to add at least one boundary.\n");
// Store the field elements
for (auto e : fldGrp->tag->elems)
@@ -185,7 +169,7 @@ void MshDeform::initialize()
uniqueSort(nodes1);
// Create the nodes pair
if (nodes0.size() != nodes1.size())
- throw std::runtime_error("MshDeform: internal boudnaries must have the same number of nodes!");
+ throw std::runtime_error("MshDeform::initialize: internal boudnaries must have the same number of nodes!");
for (auto n0 : nodes0)
{
for (auto n1 : nodes1)
@@ -301,14 +285,12 @@ Eigen::MatrixXd MshDeform::buildK(tbox::Element const &e, Eigen::MatrixXd const
*/
void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
{
- // Multithread
- tbb::spin_mutex mutex;
-
// List of triplets to build stifness matrix
std::deque<Eigen::Triplet<double>> T;
// Build stiffness matrix
- tbb::parallel_for_each(fldElems.begin(), fldElems.end(), [&](Element *e) {
+ for (Element const *e : fldElems)
+ {
// Hooke tensor
Eigen::MatrixXd H;
if (nDim == 2)
@@ -345,7 +327,6 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
// Elementary stiffness matrix
Eigen::MatrixXd Ke = this->buildK(*e, H);
// Assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
{
for (int m = 0; m < nDim; ++m)
@@ -361,7 +342,7 @@ void MshDeform::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &K)
}
}
}
- });
+ }
// Periodic boundary condition step 2
// -> for each node pair(a, b): write "pos_a - pos_b" = 0 on row a
for (auto p : intBndPair)
@@ -446,7 +427,7 @@ void MshDeform::build(Eigen::VectorXd &f)
*/
void MshDeform::deform()
{
- // Multithread
+ // Multithread (required for linear solver)
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
// Stiffness matrix
@@ -455,28 +436,24 @@ void MshDeform::deform()
Eigen::VectorXd f = Eigen::VectorXd::Zero(nDim * mshSize), q = Eigen::VectorXd::Zero(nDim * mshSize);
// Build stiffness matrix and RHS
- std::cout << "- Deforming mesh:\n"
- << "assembling matrices... " << std::flush;
+ if (verbose > 0)
+ std::cout << "- Deforming mesh:\n"
+ << "assembling matrices... " << std::flush;
this->build(K);
this->build(f);
- std::cout << "done" << std::endl;
+ if (verbose > 0)
+ std::cout << "done" << std::endl;
// Solve the SoE
- std ::cout << "solving equations... " << std::flush;
- Gmres gmres;
- // set the preconditioner
- gmres.setFillFactor(2);
- gmres.setDropTol(1e-6);
- // set the gmres
- gmres.setTolerance(1e-8);
- gmres.setRestart(50);
- gmres.setMaxIterations(500);
- // solve
+ if (verbose > 0)
+ std ::cout << "solving equations... " << std::flush;
Eigen::Map<Eigen::VectorXd> f_(f.data(), f.size()), q_(q.data(), q.size());
- gmres.compute(K, f_, q_);
- std::cout << "done (#it: " << gmres.getIterations() << ", error: " << std::scientific << gmres.getError() << std::fixed << ")" << std::endl;
+ linsol->compute(K, f_, q_);
+ if (verbose > 0)
+ std::cout << "done (#it: " << linsol->getIterations() << ", error: " << std::scientific << linsol->getError() << std::fixed << ")" << std::endl;
// Update position
- std ::cout << "updating mesh... " << std::flush;
+ if (verbose > 0)
+ std ::cout << "updating mesh... " << std::flush;
for (auto n : msh->nodes)
{
for (int m = 0; m < nDim; m++)
@@ -497,8 +474,9 @@ void MshDeform::deform()
e->update();
for (auto e : intBndElems)
e->update();
- std::cout << "done\n"
- << std::endl;
+ if (verbose > 0)
+ std::cout << "done\n"
+ << std::endl;
}
/**
diff --git a/tbox/src/wMshDeform.h b/tbox/src/wMshDeform.h
index e968058397e19a090196e78dadbe758523a59ad0..3d08c5709556b13364a09afbe17728232aff9dcc 100644
--- a/tbox/src/wMshDeform.h
+++ b/tbox/src/wMshDeform.h
@@ -36,10 +36,9 @@ class TBOX_API MshDeform : public fwk::wSharedObject
private:
size_t mshSize; ///< number of nodes in the mesh
int nDim; ///< dimension of the problem (2 or 3)
+ int nthreads; ///< number of threads
+ int verbose; ///< verbosity level
std::vector<int> rows; ///< unknown local index
- bool field; ///< flag to check if field has been added
- bool fixed; ///< flag to check if at least one fixed boundary has been added
- bool moving; ///< flag to check if at least one moving boundary has been added
Group *fldGrp; ///< temporary group for field
std::vector<Group *> fxdBndGrp; ///< temporary group for fixed boundaries
std::vector<Group *> movBndGrp; ///< temporary group for moving boundaries
@@ -60,17 +59,17 @@ private:
Eigen::MatrixXd buildK(tbox::Element const &e, Eigen::MatrixXd const &H);
public:
- std::shared_ptr<MshData> msh;
- int nthreads;
+ std::shared_ptr<MshData> msh; ///< mesh
+ std::shared_ptr<tbox::LinearSolver> linsol; ///< linear (inner) solver
- MshDeform(std::shared_ptr<MshData> _msh, int _nDim, int nthrds = 1);
+ MshDeform(std::shared_ptr<MshData> _msh, std::shared_ptr<tbox::LinearSolver> _linsol, int _nDim, int nthrds = 1, int vrb = 1);
virtual ~MshDeform() { std::cout << "~MshDeform()\n"; }
void setField(std::string const &fldName);
- void addMoving(std::vector<std::string> const &fxdBndName);
- void addFixed(std::vector<std::string> const &movBndName);
+ void addMoving(std::string const &fxdBndName);
+ void addFixed(std::string const &movBndName);
void setSymmetry(std::string const &symBndName, size_t _blck);
- void addInternal(std::vector<std::string> const &intBndName);
+ void addInternal(std::string const &intBndName, std::string const &intBndName_);
void initialize();
void savePos();
@@ -84,4 +83,4 @@ public:
} // namespace tbox
-#endif //WMSHDEFORM_H
\ No newline at end of file
+#endif // WMSHDEFORM_H
\ No newline at end of file
diff --git a/tbox/src/wMshExport.h b/tbox/src/wMshExport.h
index 50627aa1df3ce090c589b2853ced287dbaa5449d..19d6e3a056fac6a8092d599a5c1cc5715576d858 100644
--- a/tbox/src/wMshExport.h
+++ b/tbox/src/wMshExport.h
@@ -25,6 +25,17 @@
namespace tbox
{
+/**
+ * @brief Writer type
+ * @authors Adrien Crovato
+ */
+enum class WrtType
+{
+ UNDEFINED = 0,
+ GMSH = 1,
+ VTK = 2
+};
+
/**
* @brief Base class to write mesh
* @authors Adrien Crovato, Romain Boman
@@ -38,6 +49,7 @@ public:
MshExport(std::shared_ptr<MshData> _msh);
virtual ~MshExport() {}
+ virtual WrtType type() const { return WrtType::UNDEFINED; }
virtual void save(std::string const &fname) const;
#ifndef SWIG
@@ -47,4 +59,4 @@ public:
} // namespace tbox
-#endif //WMSHEXPORT_H
+#endif // WMSHEXPORT_H
diff --git a/tbox/src/wMumps.h b/tbox/src/wMumps.h
index 479d29f551fc5717078f53ff3269516d115a9191..535aa6f3a77f757b016c09a67eac15f1bdd230e3 100644
--- a/tbox/src/wMumps.h
+++ b/tbox/src/wMumps.h
@@ -53,6 +53,7 @@ private:
public:
Mumps();
~Mumps();
+ virtual LSolType type() const override { return LSolType::MUMPS; }
#ifndef SWIG
virtual void analyze(Eigen::SparseMatrix<double> const &A) override;
diff --git a/tbox/src/wPardiso.cpp b/tbox/src/wPardiso.cpp
index 01a589036b406527c7995d4d20c940fb83871dd8..dc9b8201f5213ff70e544d98c650c563a8f12def 100644
--- a/tbox/src/wPardiso.cpp
+++ b/tbox/src/wPardiso.cpp
@@ -50,7 +50,7 @@ void Pardiso::setOption(int k, int v)
void Pardiso::write(std::ostream &out) const
{
- out << "Pardiso (MKL/Eigen)" << std::endl;
+ out << "Eigen Intel MKL PARDISO" << std::endl;
}
#endif // USE_MKL
diff --git a/tbox/src/wPardiso.h b/tbox/src/wPardiso.h
index 5936fe69e72896615cf63cfddcc1185c53351cd5..3e6c57657b4cded842547860c4c336036443ba66 100644
--- a/tbox/src/wPardiso.h
+++ b/tbox/src/wPardiso.h
@@ -38,6 +38,7 @@ class TBOX_API Pardiso : public LinearSolver
public:
Pardiso() : LinearSolver() {}
virtual ~Pardiso() {}
+ virtual LSolType type() const override { return LSolType::PARDISO; }
#ifndef SWIG
virtual void analyze(Eigen::SparseMatrix<double> const &A) override;
diff --git a/tbox/src/wPoint1.h b/tbox/src/wPoint1.h
index fc78eb37711c4af8828be1100b83465ded5e8121..b905db0db099214b61ce81a55147316d8b716465 100644
--- a/tbox/src/wPoint1.h
+++ b/tbox/src/wPoint1.h
@@ -31,7 +31,7 @@ class TBOX_API Point1 : public Element
{
public:
Point1(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
- virtual ELTYPE type() const override { return ELTYPE::POINT1; }
+ virtual ElType type() const override { return ElType::POINT1; }
private:
static CachePoint1 &getCache(int n);
diff --git a/tbox/src/wQuad4.h b/tbox/src/wQuad4.h
index 75cde412625a5cca2d1fb921bf209eef3e942006..e50f1213e4b8872815ab6fc8392ff4fa6a5f6ec2 100644
--- a/tbox/src/wQuad4.h
+++ b/tbox/src/wQuad4.h
@@ -47,7 +47,7 @@ protected:
public:
Quad4(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
virtual ~Quad4() {}
- virtual ELTYPE type() const override { return ELTYPE::QUAD4; }
+ virtual ElType type() const override { return ElType::QUAD4; }
#ifndef SWIG
virtual void initValues(bool isvol, int ordr = 2) override;
diff --git a/tbox/src/wSparseLu.cpp b/tbox/src/wSparseLu.cpp
index e22a5271e90e83c8a00158b501a11a2b0ba30c3c..50bfb308e5a083ff0fd8a55c0e2905b1fd3de2e7 100644
--- a/tbox/src/wSparseLu.cpp
+++ b/tbox/src/wSparseLu.cpp
@@ -37,5 +37,5 @@ void SparseLu::compute(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::V
void SparseLu::write(std::ostream &out) const
{
- out << "SparseLU (Eigen)" << std::endl;
+ out << "Eigen SparseLU" << std::endl;
}
diff --git a/tbox/src/wSparseLu.h b/tbox/src/wSparseLu.h
index 323a5d628b90e522114ed53e33f46f1c29ae9345..84e7abef39b0d3a6cc661840f635f4cbc0f23ed3 100644
--- a/tbox/src/wSparseLu.h
+++ b/tbox/src/wSparseLu.h
@@ -38,6 +38,7 @@ class TBOX_API SparseLu : public LinearSolver
public:
SparseLu() : LinearSolver() {}
virtual ~SparseLu() {}
+ virtual LSolType type() const override { return LSolType::SPARSE_LU; }
#ifndef SWIG
virtual void analyze(Eigen::SparseMatrix<double> const &A) override;
diff --git a/tbox/src/wTag.cpp b/tbox/src/wTag.cpp
index e0089e7924db7deb0f498ecc93f7661c818fbdc7..acc815132f80b31b1f1f9a18a272448b0c7274bb 100644
--- a/tbox/src/wTag.cpp
+++ b/tbox/src/wTag.cpp
@@ -26,7 +26,7 @@ void Tag::write(std::ostream &out) const
out << "D) named \"" << name << "\" (" << elems.size() << " elements) ";
// count each type
- std::map<ELTYPE, int> types;
+ std::map<ElType, int> types;
for (auto e : elems)
types[e->type()] += 1;
diff --git a/tbox/src/wTetra4.h b/tbox/src/wTetra4.h
index 371dcba8ae2b81b546be302805081cc297c253f0..af669380e3faba104dfe7cb4aab4b133e822bce5 100644
--- a/tbox/src/wTetra4.h
+++ b/tbox/src/wTetra4.h
@@ -45,7 +45,7 @@ protected:
public:
Tetra4(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
virtual ~Tetra4() {}
- virtual ELTYPE type() const override { return ELTYPE::TETRA4; }
+ virtual ElType type() const override { return ElType::TETRA4; }
#ifndef SWIG
virtual void initValues(bool isvol, int ordr = 2) override;
diff --git a/tbox/src/wTri3.h b/tbox/src/wTri3.h
index c042c4c5a93094652fadab0e549d35c29fb4c9c3..6dab1db7f1147760c94a5a7fd000710f75b5cf24 100644
--- a/tbox/src/wTri3.h
+++ b/tbox/src/wTri3.h
@@ -51,7 +51,7 @@ protected:
public:
Tri3(int n, Tag *_ptag, Tag *_etag, std::vector<Tag *> &_parts, std::vector<Node *> &nods);
virtual ~Tri3() {}
- virtual ELTYPE type() const override { return ELTYPE::TRI3; }
+ virtual ElType type() const override { return ElType::TRI3; }
#ifndef SWIG
virtual void initValues(bool isvol, int ordr = 2) override;
diff --git a/tbox/src/wUnitTest.cpp b/tbox/src/wUnitTest.cpp
index 7b886c29e77e50893e9d5c64991e1b34e538b72d..cc91f80460c03b11a97ac21f17223b3fd3ff41df 100644
--- a/tbox/src/wUnitTest.cpp
+++ b/tbox/src/wUnitTest.cpp
@@ -82,7 +82,7 @@ void UnitTest::mshDeform(MshDeform *mshDef, int d)
{
// Initialize memory (volume element only)
for (auto e : mshDef->msh->elems)
- if ((d == 2 && e->type() == ELTYPE::TRI3) || (d == 3 && e->type() == ELTYPE::TETRA4))
+ if ((d == 2 && e->type() == ElType::TRI3) || (d == 3 && e->type() == ElType::TETRA4))
e->initValues(true);
// Call deform
mshDef->deform();
diff --git a/tbox/tests/meshDeformation.py b/tbox/tests/meshDeformation.py
index bfd2010f4832468f3b44e9bdc7b69be0a025d0bd..78f8353c8fee407ef33e4a7b55153b33366fcc60 100644
--- a/tbox/tests/meshDeformation.py
+++ b/tbox/tests/meshDeformation.py
@@ -39,10 +39,10 @@ def main():
msh = gmsh.MeshLoader("models/beam.geo",__file__).execute(**pars)
gmshWriter = tbox.GmshExport(msh)
# create mesh deformation handler
- mshDef = tbox.MshDeform(msh, 2, nthrds=parseargs().k)
+ mshDef = tbox.MshDeform(msh, tbox.SparseLu(), 2, nthrds=parseargs().k)
mshDef.setField("internalField")
- mshDef.addFixed(["clamp"])
- mshDef.addMoving(["tip"])
+ mshDef.addFixed("clamp")
+ mshDef.addMoving("tip")
mshDef.initialize()
# deform the mesh (dummy translation of "tip" nodes)
diff --git a/tbox/tests/meshDeformation3.py b/tbox/tests/meshDeformation3.py
index 7c0599d019ff8bd58e94602b6361368f66b0face..08e67b6b30f5bd0a453a2785a0a0aab61e658f28 100644
--- a/tbox/tests/meshDeformation3.py
+++ b/tbox/tests/meshDeformation3.py
@@ -38,10 +38,10 @@ def main():
msh = gmsh.MeshLoader("models/beam3.geo",__file__).execute(**pars)
gmshWriter = tbox.GmshExport(msh)
# create mesh deformation handler
- mshDef = tbox.MshDeform(msh, 3, nthrds=parseargs().k)
+ mshDef = tbox.MshDeform(msh, tbox.SparseLu(), 3, nthrds=parseargs().k)
mshDef.setField("field")
- mshDef.addFixed(["clamp"])
- mshDef.addMoving(["tip"])
+ mshDef.addFixed("clamp")
+ mshDef.addMoving("tip")
mshDef.initialize()
# deform the mesh (dummy translation of "tip" nodes)
diff --git a/tbox/utils.py b/tbox/utils.py
index 7d32f8dc690809a80540655a25e8fec2f554b97d..73ec554d506fbd70468392746b5890de2e38da2c 100644
--- a/tbox/utils.py
+++ b/tbox/utils.py
@@ -2,13 +2,13 @@
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
-#
+#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
-#
+#
# http://www.apache.org/licenses/LICENSE-2.0
-#
+#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -23,7 +23,6 @@ Adrien Crovato
def myrange(start, stop, step):
"""Extended range with custom step
- Adrien Crovato
"""
vec = []
i = start
@@ -32,33 +31,34 @@ def myrange(start, stop, step):
i += step
return vec
-def plot(x, y, xlbl, ylbl, title, invert=False):
+def plot(x, y, cfg):
"""Plot using matplotlib
- Adrien Crovato
- @todo extend to pass unlimited number of arguments
"""
from fwk.wutils import parseargs
args = parseargs()
if not args.nogui:
import matplotlib.pyplot as plt
# define font
- font = {'family': 'serif',
- 'color': 'darkred',
- 'weight': 'normal',
- 'size': 16,}
+ font = {'family': 'serif', 'color': 'darkred', 'weight': 'normal', 'size': 16}
# plot results
- plt.plot(x, y, 'b-',lw=2)
- plt.xlabel(xlbl, fontdict=font)
- plt.ylabel(ylbl, fontdict=font)
- plt.title(title, fontdict=font)
- if invert:
+ plt.plot(x, y, 'b-', lw=2)
+ if 'xlabel' in cfg:
+ plt.xlabel(cfg['xlabel'], fontdict=font)
+ if 'ylabel' in cfg:
+ plt.ylabel(cfg['ylabel'], fontdict=font)
+ if 'title' in cfg:
+ plt.title(cfg['title'], fontdict=font)
+ if 'xlim' in cfg:
+ plt.xlim(cfg['xlim'])
+ if 'ylim' in cfg:
+ plt.ylim(cfg['ylim'])
+ if 'invert' in cfg and cfg['invert']:
plt.gca().invert_yaxis()
- #plt.grid(True)
+ # display
plt.show()
def sort(id, data):
"""Sort data array against id (line connectivity list)
- Adrien Crovato
"""
import numpy as np
# sort id vector
@@ -76,7 +76,6 @@ def sort(id, data):
def read(filename):
"""Read from file and store in data array
- Adrien Croavto
"""
import io
import numpy as np
@@ -89,7 +88,6 @@ def read(filename):
def write(data, name, frmt, dlm, hdr, cmts):
"""Write data array to file
- Adrien Crovato
"""
import numpy as np
print('writing data file ' + name + '.')
diff --git a/tbox/viewer.py b/tbox/viewer.py
index ddabbc90fe01ec8dbc3c357e47f3d141c87dc3c1..c5f42697cfbbd50caac3dac6327197f64fb66acb 100644
--- a/tbox/viewer.py
+++ b/tbox/viewer.py
@@ -2,13 +2,13 @@
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
-#
+#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
-#
+#
# http://www.apache.org/licenses/LICENSE-2.0
-#
+#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
@@ -35,10 +35,13 @@ class GUI:
if not os.path.isfile(optfile):
self.defopt(optfile)
cmd="gmsh %s %s %s" % (mshfile, posfile, optfile)
- os.system(cmd)
+ os.system(cmd)
def defopt(self,fname='defaultflowOpt.opt'):
f = open(fname,'w')
f.write("General.Orthographic = 0;\n")
+ f.write('View[0].IntervalsType = 3;\n')
+ f.write('For i In {1:PostProcessing.NbViews-1}\n')
+ f.write(' View[i].Visible = 0;\n')
+ f.write('EndFor\n')
f.close()
-
diff --git a/tboxVtk/_src/CMakeLists.txt b/tboxVtk/_src/CMakeLists.txt
index f1ad6f73a27613af52c9e94bbd1f84d55285a4f0..0011635748f147cd815dc45199cc4ad1a955ad05 100644
--- a/tboxVtk/_src/CMakeLists.txt
+++ b/tboxVtk/_src/CMakeLists.txt
@@ -19,35 +19,18 @@ INCLUDE(${SWIG_USE_FILE})
FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
FILE(GLOB ISRCS *.i)
-SET(CMAKE_SWIG_FLAGS "")
SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} "-interface" "_tboxVtkw") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(tboxVtkw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET tboxVtkw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+MACRO_DebugPostfix(tboxVtkw)
-SET(SWINCFLAGS
--I${PROJECT_SOURCE_DIR}/tboxVtk/src
--I${PROJECT_SOURCE_DIR}/tbox/src
--I${PROJECT_SOURCE_DIR}/tbox/_src
--I${PROJECT_SOURCE_DIR}/fwk/src
--I${PROJECT_SOURCE_DIR}/fwk/_src
-)
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
-
-if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
- SWIG_ADD_MODULE(tboxVtkw python ${ISRCS} ${SRCS})
-else()
- SWIG_ADD_LIBRARY(tboxVtkw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-endif()
-MACRO_DebugPostfix(_tboxVtkw)
-
-TARGET_INCLUDE_DIRECTORIES(_tboxVtkw PRIVATE ${PROJECT_SOURCE_DIR}/tboxVtk/_src
- ${PROJECT_SOURCE_DIR}/fwk/_src
- ${PROJECT_SOURCE_DIR}/tbox/_src
- ${PYTHON_INCLUDE_PATH}
-)
-
-SWIG_LINK_LIBRARIES(tboxVtkw
- tboxVtk tbox fwk ${PYTHON_LIBRARIES}
-)
+TARGET_INCLUDE_DIRECTORIES(tboxVtkw PRIVATE ${PROJECT_SOURCE_DIR}/tboxVtk/_src
+ ${PROJECT_SOURCE_DIR}/tbox/_src
+ ${PROJECT_SOURCE_DIR}/fwk/_src
+ ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(tboxVtkw PRIVATE tboxVtk tbox fwk ${PYTHON_LIBRARIES})
INSTALL(FILES "${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/tboxVtkw.py"
DESTINATION ${CMAKE_INSTALL_PREFIX})
-INSTALL(TARGETS _tboxVtkw DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(TARGETS tboxVtkw DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/tboxVtk/src/CMakeLists.txt b/tboxVtk/src/CMakeLists.txt
index 2d38e5c6013bd1b04b3e5e32b1f2cdbb03524005..288a3a4a674d94056a5c2ae8924bbbe1047596d7 100644
--- a/tboxVtk/src/CMakeLists.txt
+++ b/tboxVtk/src/CMakeLists.txt
@@ -27,14 +27,18 @@ TARGET_INCLUDE_DIRECTORIES(tboxVtk PUBLIC ${PROJECT_SOURCE_DIR}/tboxVtk/src)
# 1. first, look for all available VTK libraries
-FIND_PACKAGE(VTK REQUIRED NO_MODULE)
+FIND_PACKAGE(VTK REQUIRED NO_MODULE QUIET)
IF(${VTK_VERSION} VERSION_LESS 6.0)
MESSAGE(FATAL_ERROR "VTK>=6 required!")
ENDIF()
-MESSAGE(STATUS "VTK_LIBRARIES1 = ${VTK_LIBRARIES}") # display all VTK libs
+#MESSAGE(STATUS "VTK_LIBRARIES1 = ${VTK_LIBRARIES}") # display all VTK libs
# list of required components
-SET(REQVTKCMP vtkIOXML)
+IF(${VTK_VERSION} VERSION_LESS 9.0)
+ SET(REQVTKCMP vtkIOXML)
+ELSE()
+ SET(REQVTKCMP IOXML)
+ENDIF()
# check whether MPI has been compiled
IF(TARGET vtkParallelMPI) # NOTE: this test also works with a previous single-COMPONENT search
@@ -47,16 +51,19 @@ ENDIF()
# 2. second, perform a new search with a limited number of components
# NOTE: you may just disable the following line to link against ALL vtk libraries!
-FIND_PACKAGE(VTK REQUIRED COMPONENTS ${REQVTKCMP} NO_MODULE)
+FIND_PACKAGE(VTK REQUIRED COMPONENTS ${REQVTKCMP} NO_MODULE QUIET)
+IF(${VTK_FOUND})
+ MESSAGE(STATUS "Found VTK ${VTK_VERSION} at \"${VTK_PREFIX_PATH}\" (components: ${VTK_LIBRARIES})")
+ENDIF()
# debug
-MESSAGE(STATUS "VTK_DEFINITIONS = ${VTK_DEFINITIONS}")
-MESSAGE(STATUS "VTK_REQUIRED_CXX_FLAGS = ${VTK_REQUIRED_CXX_FLAGS}")
-MESSAGE(STATUS "VTK_REQUIRED_EXE_LINKER_FLAGS = ${VTK_REQUIRED_EXE_LINKER_FLAGS}")
-MESSAGE(STATUS "VTK_REQUIRED_SHARED_LINKER_FLAGS = ${VTK_REQUIRED_SHARED_LINKER_FLAGS}")
-MESSAGE(STATUS "VTK_USE_FILE = ${VTK_USE_FILE}")
-MESSAGE(STATUS "VTK_INCLUDE_DIR = ${VTK_INCLUDE_DIR}")
-MESSAGE(STATUS "VTK_LIBRARIES = ${VTK_LIBRARIES}")
+#MESSAGE(STATUS "VTK_DEFINITIONS = ${VTK_DEFINITIONS}")
+#MESSAGE(STATUS "VTK_REQUIRED_CXX_FLAGS = ${VTK_REQUIRED_CXX_FLAGS}")
+#MESSAGE(STATUS "VTK_REQUIRED_EXE_LINKER_FLAGS = ${VTK_REQUIRED_EXE_LINKER_FLAGS}")
+#MESSAGE(STATUS "VTK_REQUIRED_SHARED_LINKER_FLAGS = ${VTK_REQUIRED_SHARED_LINKER_FLAGS}")
+#MESSAGE(STATUS "VTK_USE_FILE = ${VTK_USE_FILE}")
+#MESSAGE(STATUS "VTK_INCLUDE_DIR = ${VTK_INCLUDE_DIR}")
+#MESSAGE(STATUS "VTK_LIBRARIES = ${VTK_LIBRARIES}")
TARGET_INCLUDE_DIRECTORIES(tboxVtk PUBLIC ${VTK_INCLUDE_DIRS})
TARGET_COMPILE_DEFINITIONS(tboxVtk PUBLIC ${VTK_DEFINITIONS})
diff --git a/tboxVtk/src/wVtkExport.cpp b/tboxVtk/src/wVtkExport.cpp
index 322eff828e5a4f08080f0a15fecf2c83e8caad22..26d1750fae91e2d73c5059d475dca71a092f5d92 100644
--- a/tboxVtk/src/wVtkExport.cpp
+++ b/tboxVtk/src/wVtkExport.cpp
@@ -190,42 +190,42 @@ void VtkExport::create()
// Create elements
for (auto e : msh->elems)
{
- if (e->type() == ELTYPE::HEX8)
+ if (e->type() == ElType::HEX8)
{
vtkHexaP cell = vtkHexaP::New();
for (size_t i = 0; i < e->nodes.size(); ++i)
cell->GetPointIds()->SetId(i, e->nodes[i]->no - 1);
grid->InsertNextCell(cell->GetCellType(), cell->GetPointIds());
}
- else if (e->type() == ELTYPE::TETRA4)
+ else if (e->type() == ElType::TETRA4)
{
vtkTetraP cell = vtkTetraP::New();
for (size_t i = 0; i < e->nodes.size(); ++i)
cell->GetPointIds()->SetId(i, e->nodes[i]->no - 1);
grid->InsertNextCell(cell->GetCellType(), cell->GetPointIds());
}
- else if (e->type() == ELTYPE::QUAD4)
+ else if (e->type() == ElType::QUAD4)
{
vtkQuadP cell = vtkQuadP::New();
for (size_t i = 0; i < e->nodes.size(); ++i)
cell->GetPointIds()->SetId(i, e->nodes[i]->no - 1);
grid->InsertNextCell(cell->GetCellType(), cell->GetPointIds());
}
- else if (e->type() == ELTYPE::TRI3)
+ else if (e->type() == ElType::TRI3)
{
vtkTriP cell = vtkTriP::New();
for (size_t i = 0; i < e->nodes.size(); ++i)
cell->GetPointIds()->SetId(i, e->nodes[i]->no - 1);
grid->InsertNextCell(cell->GetCellType(), cell->GetPointIds());
}
- else if (e->type() == ELTYPE::LINE2)
+ else if (e->type() == ElType::LINE2)
{
vtkLineP cell = vtkLineP::New();
for (size_t i = 0; i < e->nodes.size(); ++i)
cell->GetPointIds()->SetId(i, e->nodes[i]->no - 1);
grid->InsertNextCell(cell->GetCellType(), cell->GetPointIds());
}
- else if (e->type() == ELTYPE::POINT1)
+ else if (e->type() == ElType::POINT1)
{
vtkVtxP cell = vtkVtxP::New();
for (size_t i = 0; i < e->nodes.size(); ++i)
diff --git a/tboxVtk/src/wVtkExport.h b/tboxVtk/src/wVtkExport.h
index 8dea4d0fdfd78c8a35c0ee1db8991a4e7ce4932b..eca8b73cb475956294e8631e6621c12d852c1cda 100644
--- a/tboxVtk/src/wVtkExport.h
+++ b/tboxVtk/src/wVtkExport.h
@@ -44,6 +44,7 @@ class TBOXVTK_API VtkExport : public tbox::MshExport
public:
VtkExport(std::shared_ptr<tbox::MshData> _msh);
virtual ~VtkExport();
+ virtual tbox::WrtType type() const override { return tbox::WrtType::VTK; }
virtual void save(std::string const &fname) const override;
#ifndef SWIG
@@ -55,4 +56,4 @@ public:
} // namespace tboxVtk
-#endif //WGMSHUTILS_H
+#endif // WVTKEXPORT_H