diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000000000000000000000000000000000000..60cf9d536e56d8675468717814ffb5698298cf1e
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,11 @@
+---
+# https://clang.llvm.org/docs/ClangFormatStyleOptions.html#
+# the "Visual Studio" style is similar to:
+BasedOnStyle: LLVM
+UseTab: Never
+IndentWidth: 4
+BreakBeforeBraces: Allman
+AccessModifierOffset: -4
+SortIncludes: false
+ColumnLimit: 0
+...
diff --git a/.editorconfig b/.editorconfig
new file mode 100644
index 0000000000000000000000000000000000000000..de6044b96d078e062de7f658aef508ff77461564
--- /dev/null
+++ b/.editorconfig
@@ -0,0 +1,8 @@
+# EditorConfig is awesome: https://EditorConfig.org
+
+root = true
+
+[*]
+charset = utf-8
+indent_style = space
+indent_size = 4
diff --git a/.gitignore b/.gitignore
index 109a5c8905a52423e8e2214a2c8b8621554f8cf6..4bb53200ef03cc1a3c0738d0797d1a5a2034923e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,4 @@ __pycache__/
*.zip
*.7z
*.rar
+build/
diff --git a/cxxfem/CMakeLists.txt b/cxxfem/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..be50e6426998b5bab4d01a37e76b8461ff3965db
--- /dev/null
+++ b/cxxfem/CMakeLists.txt
@@ -0,0 +1,128 @@
+PROJECT(FEM CXX)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.14)
+
+OPTION(FEM_USE_MKL "Use Pardiso linear solver from MKL" ON)
+
+
+# ------------------------------------------------------------------------------
+# Find libraries and setup compiler
+# ------------------------------------------------------------------------------
+
+# put executables/libraries together in bin/ folder
+SET(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin CACHE PATH
+ "Single output directory for building all libraries.")
+SET(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin CACHE PATH
+ "Single output directory for building all executables.")
+MARK_AS_ADVANCED(CMAKE_LIBRARY_OUTPUT_DIRECTORY
+ CMAKE_RUNTIME_OUTPUT_DIRECTORY)
+
+# build type is "" by default in Linux
+IF(NOT CMAKE_BUILD_TYPE)
+ SET( CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE)
+ENDIF()
+
+IF(CMAKE_CXX_COMPILER_ID MATCHES "MSVC")
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /MP") # parallel build with MSVC
+ENDIF()
+
+# enable C++11
+SET(CMAKE_CXX_STANDARD 17)
+SET(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+IF(APPLE)
+ # on macOS, do not give priority to frameworks/apps
+ SET(CMAKE_FIND_APPBUNDLE LAST)
+ SET(CMAKE_FIND_FRAMEWORK LAST)
+ENDIF()
+
+# find gmsh-sdk
+# gmsh.h
+FIND_PATH(GMSH_INCLUDE_DIRS NAMES "gmsh.h")
+MESSAGE(STATUS "GMSH_INCLUDE_DIRS=" ${GMSH_INCLUDE_DIRS})
+if(NOT GMSH_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "gmsh.h not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${GMSH_INCLUDE_DIRS})
+
+# libgmsh.so
+FIND_LIBRARY(GMSH_LIBRARIES gmsh)
+MESSAGE(STATUS "GMSH_LIBRARIES=" ${GMSH_LIBRARIES})
+IF(NOT GMSH_LIBRARIES)
+ MESSAGE(FATAL_ERROR "gmsh library not found!")
+ENDIF()
+
+# find MKL
+
+IF(FEM_USE_MKL)
+ # header (mkl.h) searched using INCLUDE or MKLROOT
+ FIND_PATH(MKL_INCLUDE_DIRS NAMES "mkl.h" PATHS "$ENV{MKLROOT}/include")
+ MESSAGE(STATUS "MKL_INCLUDE_DIRS=${MKL_INCLUDE_DIRS}")
+ # library (mkl_rt.so) searched using LIBRARY_PATH (Linux/macOS) or LIB (windows)
+ FIND_LIBRARY(MKL_LIBRARIES mkl_rt PATHS ENV LIBRARY_PATH)
+ IF(MINGW)
+ FIND_LIBRARY(MKL_LIB1 mkl_intel_lp64_dll PATHS ENV LIBRARY_PATH)
+ FIND_LIBRARY(MKL_LIB2 mkl_intel_thread_dll PATHS ENV LIBRARY_PATH)
+ FIND_LIBRARY(MKL_LIB3 mkl_core_dll PATHS ENV LIBRARY_PATH)
+ FIND_LIBRARY(MKL_LIB4 libiomp5md PATHS ENV LIBRARY_PATH)
+ # MESSAGE(STATUS "MKL_LIB1=${MKL_LIB1}")
+ # MESSAGE(STATUS "MKL_LIB2=${MKL_LIB2}")
+ # MESSAGE(STATUS "MKL_LIB3=${MKL_LIB3}")
+ # MESSAGE(STATUS "MKL_LIB4=${MKL_LIB4}")
+ SET(MKL_LIBRARIES ${MKL_LIB1} ${MKL_LIB2} ${MKL_LIB3} ${MKL_LIB4})
+ ENDIF()
+ MESSAGE(STATUS "MKL_LIBRARIES=${MKL_LIBRARIES}")
+ENDIF()
+
+# find Eigen
+FIND_PATH(EIGEN_INCLUDE_DIRS "Eigen/Dense"
+ PATHS "${PROJECT_SOURCE_DIR}/lib/eigen" "/usr/include/eigen3")
+MESSAGE(STATUS "EIGEN_INCLUDE_DIRS=" ${EIGEN_INCLUDE_DIRS})
+IF(NOT EIGEN_INCLUDE_DIRS)
+ MESSAGE(FATAL_ERROR "Eigen include dir not found!")
+ENDIF()
+INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS})
+
+# increase warning level
+INCLUDE(CheckCXXCompilerFlag)
+CHECK_CXX_COMPILER_FLAG("-Wall" WALL)
+IF(WALL AND NOT MSVC)
+ ADD_COMPILE_OPTIONS(-Wall)
+ENDIF()
+
+# OpenMP --
+FIND_PACKAGE(OpenMP)
+IF(OPENMP_FOUND)
+ MESSAGE ("OpenMP enabled (CXX_FLAGS=${OpenMP_CXX_FLAGS})")
+ELSE()
+ MESSAGE ("OpenMP disabled.")
+ENDIF()
+
+# Python --
+# use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
+FIND_PACKAGE(Python3 COMPONENTS Interpreter Development REQUIRED)
+MESSAGE(STATUS "Python3_EXECUTABLE=${Python3_EXECUTABLE}")
+MESSAGE(STATUS "Python3_LIBRARIES=${Python3_LIBRARIES}")
+MESSAGE(STATUS "Python3_INCLUDE_DIRS=${Python3_INCLUDE_DIRS}")
+MESSAGE(STATUS "Python3_VERSION=${Python3_VERSION}")
+
+# SWIG --
+FIND_PACKAGE(SWIG REQUIRED)
+IF(CMAKE_GENERATOR MATCHES "Visual Studio") # not MSVC because of nmake & jom
+ SET(CMAKE_SWIG_OUTDIR "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/$(Configuration)/")
+ELSE()
+ SET(CMAKE_SWIG_OUTDIR "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}")
+ENDIF()
+SET(CMAKE_SWIG_FLAGS -py3 -O)
+INCLUDE(${SWIG_USE_FILE})
+
+
+
+
+# ------------------------------------------------------------------------------
+CONFIGURE_FILE(config.h.in config.h)
+INCLUDE_DIRECTORIES(${PROJECT_BINARY_DIR})
+
+ADD_SUBDIRECTORY( src )
+ADD_SUBDIRECTORY( exe )
+ADD_SUBDIRECTORY( _src )
+
diff --git a/cxxfem/__init__.py b/cxxfem/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a3577e00eb47acd97f625cc658016a766d30094
--- /dev/null
+++ b/cxxfem/__init__.py
@@ -0,0 +1,4 @@
+# -*- coding: utf-8 -*-
+# fem MODULE initialization file
+
+from femi import *
diff --git a/cxxfem/_src/CMakeLists.txt b/cxxfem/_src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f3549cc40fadd9fdb29ddb5bf1dda8023b028f5f
--- /dev/null
+++ b/cxxfem/_src/CMakeLists.txt
@@ -0,0 +1,17 @@
+# CMake input file of the SWIG wrapper around "fem.so"
+# Uses new policy CMP0078 (UseSWIG generates standard target names)
+# Requires CMake >=3.13
+
+FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
+FILE(GLOB ISRCS *.i)
+
+SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS "-interface" "_femi") # avoids "import _module_d" with MSVC/Debug
+SWIG_ADD_LIBRARY(femi LANGUAGE python SOURCES ${ISRCS} ${SRCS})
+SET_PROPERTY(TARGET femi PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
+# MACRO_DebugPostfix(femi)
+
+TARGET_INCLUDE_DIRECTORIES(femi PRIVATE
+ ${PROJECT_SOURCE_DIR}/src
+ ${Python3_INCLUDE_DIRS})
+TARGET_LINK_LIBRARIES(femi PRIVATE fem ${Python3_LIBRARIES})
diff --git a/cxxfem/_src/femi.i b/cxxfem/_src/femi.i
new file mode 100644
index 0000000000000000000000000000000000000000..41956c6d8aee6d5576b4b3205edbfaad6e8b17b0
--- /dev/null
+++ b/cxxfem/_src/femi.i
@@ -0,0 +1,82 @@
+// SWIG input file of the 'fem' module
+
+%feature("autodoc","1");
+
+%module(docstring=
+"'femi' module: FEM
+(c) ULiège - A&M",
+threads="1"
+) femi
+%{
+
+#include "fem.h"
+#include "femDirichlet.h"
+#include "femNeumann.h"
+#include "femProblem.h"
+#include "femSolver.h"
+#include "femPost.h"
+#include "femMedium.h"
+#include "femExtractor.h"
+
+#include <string>
+// #include <sstream>
+
+%}
+
+%include "std_string.i"
+%include "std_vector.i"
+%include "std_map.i"
+
+// Instantiate templates
+namespace std {
+ %template(StringDoubleMap) std::map<std::string, double>;
+ %template(StringVector) std::vector<std::string>;
+ %template(DoubleVector) std::vector<double>;
+}
+
+// --------- EXCEPTIONS ---------------
+
+%include "exception.i"
+
+// from: http://swig.10945.n7.nabble.com/Trapping-Swig-DirectorException-td6013.html
+// le code suivant permet de voir la call stack dans les appels C++ => python
+
+%{
+ static void handle_exception(void) {
+ try {
+ throw;
+ } catch (std::exception &e) {
+ std::stringstream txt;
+ txt << e.what(); // << ' ' << typeid(e).name();
+ PyErr_SetString(PyExc_RuntimeError, e.what());
+ }
+ catch(...)
+ {
+ PyErr_SetString(PyExc_RuntimeError, "Unknown C++ Runtime Error");
+ }
+ }
+%}
+
+%exception {
+ try {
+ $action
+ } catch (...) {
+ // Note that if a director method failed, the Python error indicator
+ // already contains full details of the exception, and it will be
+ // reraised when we go to SWIG_fail; so no need to convert the C++
+ // exception back to a Python one
+ if (!PyErr_Occurred()) {
+ handle_exception();
+ }
+ SWIG_fail;
+ }
+}
+
+// ----------- FEM CLASSES ---------------
+%include "fem.h"
+%include "femDirichlet.h"
+%include "femNeumann.h"
+%include "femSolver.h"
+%include "femProblem.h"
+%include "femPost.h"
+%include "femMedium.h"
diff --git a/cxxfem/config.h.in b/cxxfem/config.h.in
new file mode 100644
index 0000000000000000000000000000000000000000..715648e7ecd248f2e16e6e270a85cd48e9530a8e
--- /dev/null
+++ b/cxxfem/config.h.in
@@ -0,0 +1,2 @@
+#cmakedefine PROJECT_SOURCE_DIR "@PROJECT_SOURCE_DIR@"
+#cmakedefine FEM_USE_MKL
diff --git a/cxxfem/exe/CMakeLists.txt b/cxxfem/exe/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e02d9e4445ceb4f395597b9f89c61271a919d6fd
--- /dev/null
+++ b/cxxfem/exe/CMakeLists.txt
@@ -0,0 +1,9 @@
+
+# compile each .cpp into an exe
+FILE(GLOB SRCS *.cpp)
+FOREACH(SFILE ${SRCS})
+ GET_FILENAME_COMPONENT(PROG_NAME ${SFILE} NAME_WE)
+ ADD_EXECUTABLE(${PROG_NAME} ${SFILE})
+ TARGET_LINK_LIBRARIES(${PROG_NAME} PUBLIC fem)
+ TARGET_INCLUDE_DIRECTORIES(fem PUBLIC ${PROJECT_SOURCE_DIR}/exe)
+ENDFOREACH()
diff --git a/cxxfem/exe/beam.cpp b/cxxfem/exe/beam.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..70bab05cc7708dfbdf936bcc8bf5d4fcf629454d
--- /dev/null
+++ b/cxxfem/exe/beam.cpp
@@ -0,0 +1,66 @@
+// build & run:
+// win:
+// set OMP_NUM_THREADS=10
+// cmake .. && make -j 10 && bin\square.exe
+// win debug:
+// cmake -DCMAKE_BUILD_TYPE=Debug .. && make -j 10 && bin\square.exe
+//
+// linux:
+// export OMP_NUM_THREADS=10
+// cmake .. && make -j 10 && ./bin/square
+// linux debug:
+// cmake -DCMAKE_BUILD_TYPE=Debug .. && make -j 10 && ./bin/square
+
+#include "fem.h"
+#include "femProblem.h"
+#include "femMedium.h"
+#include "femDirichlet.h"
+#include "femSolver.h"
+#include "femPost.h"
+#include "OpenMP.h"
+#include "build_square.h"
+#include "femGmsh.h"
+#include <stdexcept>
+
+int main(int argc, char **argv)
+{
+ try
+ {
+ // If compiled with OpenMP support (ubuntu), gmsh::initialize
+ // also sets the number of threads to "General.NumThreads".
+ int nthreads = OpenMP::get_max_threads();
+ std::cout << "nthreads = " << nthreads << '\n';
+ gmsh::initialize();
+ OpenMP::set_num_threads(nthreads);
+ gmsh::option::setNumber("General.Terminal", 1);
+
+ // build geometry & mesh
+ build_square(0.0, 0.0, 10.0, 1.0, 80, 8, true, true);
+
+ fem::Problem pbl;
+
+ new fem::Medium(pbl, "interior", 100.0, 0.3);
+
+ new fem::Dirichlet(pbl, "left", "X", 0.0);
+ new fem::Dirichlet(pbl, "left", "Y", 0.0);
+ new fem::Dirichlet(pbl, "pt3", "Y", -0.1);
+
+ fem::Solver solver(pbl);
+ solver.solve();
+
+ fem::Post post(solver);
+ post.build_views();
+ post.show_views( { "Y" });
+ post.view();
+
+ // gmsh::fltk::run();
+
+ gmsh::finalize();
+ return EXIT_SUCCESS;
+ }
+ catch (const std::exception &e)
+ {
+ std::cerr << "\n**ERROR: " << e.what() << '\n';
+ return EXIT_FAILURE;
+ }
+}
diff --git a/cxxfem/exe/build_square.h b/cxxfem/exe/build_square.h
new file mode 100644
index 0000000000000000000000000000000000000000..a5ed152276f03227624a9c793e9dc905bb06207c
--- /dev/null
+++ b/cxxfem/exe/build_square.h
@@ -0,0 +1,73 @@
+#ifndef BUILD_SQUARE_H
+#define BUILD_SQUARE_H
+
+/// build a square using Gmsh API
+
+#include "femGmsh.h"
+
+inline void build_square(double cx = 0.0, double cy = 0.0,
+ double Lx = 1.0, double Ly = 1.0,
+ int nx = 2, int ny = 2,
+ bool structured = true, bool recombine = true)
+{
+ gmsh::model::add("square");
+ double lc = 1.0; // no impact
+ gmsh::model::geo::addPoint(cx, cy, 0, lc, 1);
+ gmsh::model::geo::addPoint(cx + Lx, cy, 0, lc, 2);
+ gmsh::model::geo::addPoint(cx + Lx, cy + Ly, 0, lc, 3);
+ gmsh::model::geo::addPoint(cx, cy + Ly, 0, lc, 4);
+
+ gmsh::model::geo::addLine(1, 2, 1);
+ gmsh::model::geo::addLine(2, 3, 2);
+ gmsh::model::geo::addLine(4, 3, 3);
+ gmsh::model::geo::addLine(1, 4, 4);
+
+ gmsh::model::geo::addCurveLoop({1, 2, -3, -4}, 1);
+
+ gmsh::model::geo::addPlaneSurface({1}, 1);
+
+ gmsh::model::geo::synchronize(); // sinon warning "unknown entity..." lors du addPhysicalGroup
+
+ // mesh lines
+ gmsh::model::mesh::setTransfiniteCurve(1, nx + 1);
+ gmsh::model::mesh::setTransfiniteCurve(3, nx + 1);
+ gmsh::model::mesh::setTransfiniteCurve(2, ny + 1);
+ gmsh::model::mesh::setTransfiniteCurve(4, ny + 1);
+
+ // recombine
+ if (structured)
+ gmsh::model::mesh::setTransfiniteSurface(1);
+ if (recombine)
+ gmsh::model::mesh::setRecombine(2, 1); // quads
+
+ // physical groups - corners
+ auto pt1 = gmsh::model::addPhysicalGroup(0, {1});
+ auto pt2 = gmsh::model::addPhysicalGroup(0, {2});
+ auto pt3 = gmsh::model::addPhysicalGroup(0, {3});
+ auto pt4 = gmsh::model::addPhysicalGroup(0, {4});
+ gmsh::model::setPhysicalName(0, pt1, "pt1");
+ gmsh::model::setPhysicalName(0, pt2, "pt2");
+ gmsh::model::setPhysicalName(0, pt3, "pt3");
+ gmsh::model::setPhysicalName(0, pt4, "pt4");
+
+ // physical groups - lines
+ auto bottom = gmsh::model::addPhysicalGroup(1, {1}, 1000);
+ auto right = gmsh::model::addPhysicalGroup(1, {2});
+ auto top = gmsh::model::addPhysicalGroup(1, {3});
+ auto left = gmsh::model::addPhysicalGroup(1, {4});
+ gmsh::model::setPhysicalName(1, bottom, "bottom");
+ gmsh::model::setPhysicalName(1, right, "right");
+ gmsh::model::setPhysicalName(1, top, "top");
+ gmsh::model::setPhysicalName(1, left, "left");
+
+ // physical groups - surfaces
+ auto interior = gmsh::model::addPhysicalGroup(2, {1});
+ gmsh::model::setPhysicalName(2, interior, "interior");
+
+ // finalize and mesh the geometry
+ gmsh::model::mesh::generate(2);
+
+ // gmsh::write("square.msh");
+}
+
+#endif //BUILD_SQUARE_H
diff --git a/cxxfem/exe/simple_shear.cpp b/cxxfem/exe/simple_shear.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a2e3946a0f54adb885630cf415d45d6b006ecca
--- /dev/null
+++ b/cxxfem/exe/simple_shear.cpp
@@ -0,0 +1,66 @@
+// build & run:
+// win:
+// set OMP_NUM_THREADS=10
+// cmake .. && make -j 10 && bin\square.exe
+// win debug:
+// cmake -DCMAKE_BUILD_TYPE=Debug .. && make -j 10 && bin\square.exe
+//
+// linux:
+// export OMP_NUM_THREADS=10
+// cmake .. && make -j 10 && ./bin/square
+// linux debug:
+// cmake -DCMAKE_BUILD_TYPE=Debug .. && make -j 10 && ./bin/square
+
+#include "fem.h"
+#include "femProblem.h"
+#include "femMedium.h"
+#include "femDirichlet.h"
+#include "femSolver.h"
+#include "femPost.h"
+#include "OpenMP.h"
+#include "build_square.h"
+#include "femGmsh.h"
+#include <stdexcept>
+
+int main(int argc, char **argv)
+{
+ try
+ {
+ // If compiled with OpenMP support (ubuntu), gmsh::initialize
+ // also sets the number of threads to "General.NumThreads".
+ int nthreads = OpenMP::get_max_threads();
+ std::cout << "nthreads = " << nthreads << '\n';
+ gmsh::initialize();
+ OpenMP::set_num_threads(nthreads);
+ gmsh::option::setNumber("General.Terminal", 1);
+
+ // build geometry & mesh
+ build_square(0.0, 0.0, 1.0, 1.0, 10, 10, true, true); // structured
+
+ fem::Problem pbl;
+
+ new fem::Medium(pbl, "interior", 100., 0.3);
+
+ new fem::Dirichlet(pbl, "top", "X", 0.1);
+ new fem::Dirichlet(pbl, "top", "Y", 0.0);
+ new fem::Dirichlet(pbl, "bottom", "X", 0.0);
+ new fem::Dirichlet(pbl, "bottom", "Y", 0.0);
+
+ fem::Solver solver(pbl);
+ solver.solve();
+
+ fem::Post post(solver);
+ post.build_views();
+ post.view();
+
+ // gmsh::fltk::run();
+
+ gmsh::finalize();
+ return EXIT_SUCCESS;
+ }
+ catch (const std::exception &e)
+ {
+ std::cerr << "\n**ERROR: " << e.what() << '\n';
+ return EXIT_FAILURE;
+ }
+}
diff --git a/cxxfem/exe/square.cpp b/cxxfem/exe/square.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dae56bfb9eb8cae84a5ab5659b92a72ff09b1062
--- /dev/null
+++ b/cxxfem/exe/square.cpp
@@ -0,0 +1,65 @@
+// build & run:
+// win:
+// set OMP_NUM_THREADS=10
+// cmake .. && make -j 10 && bin\square.exe
+// win debug:
+// cmake -DCMAKE_BUILD_TYPE=Debug .. && make -j 10 && bin\square.exe
+//
+// linux:
+// export OMP_NUM_THREADS=10
+// cmake .. && make -j 10 && ./bin/square
+// linux debug:
+// cmake -DCMAKE_BUILD_TYPE=Debug .. && make -j 10 && ./bin/square
+
+#include "fem.h"
+#include "femProblem.h"
+#include "femMedium.h"
+#include "femDirichlet.h"
+#include "femSolver.h"
+#include "femPost.h"
+#include "OpenMP.h"
+#include "build_square.h"
+#include "femGmsh.h"
+#include <stdexcept>
+
+int main(int argc, char **argv)
+{
+ try
+ {
+ // If compiled with OpenMP support (ubuntu), gmsh::initialize
+ // also sets the number of threads to "General.NumThreads".
+ int nthreads = OpenMP::get_max_threads();
+ std::cout << "nthreads = " << nthreads << '\n';
+ gmsh::initialize();
+ OpenMP::set_num_threads(nthreads);
+ gmsh::option::setNumber("General.Terminal", 1);
+
+ // build geometry & mesh
+ build_square(0.0, 0.0, 1.0, 1.0, 10, 10, true, true); // structured
+
+ fem::Problem pbl;
+
+ new fem::Medium(pbl, "interior", 100.0, 0.3);
+
+ new fem::Dirichlet(pbl, "bottom", "Y", 0.0);
+ new fem::Dirichlet(pbl, "left", "X", 0.0);
+ new fem::Dirichlet(pbl, "right", "X", 1.0); // tensile prescribed displ
+
+ fem::Solver solver(pbl);
+ solver.solve();
+
+ fem::Post post(solver);
+ post.build_views();
+ post.view();
+
+ // gmsh::fltk::run();
+
+ gmsh::finalize();
+ return EXIT_SUCCESS;
+ }
+ catch (const std::exception &e)
+ {
+ std::cerr << "\n**ERROR: " << e.what() << '\n';
+ return EXIT_FAILURE;
+ }
+}
diff --git a/cxxfem/src/CMakeLists.txt b/cxxfem/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..cb63466ab1083dea6ad58cd6c302f730d5db3f35
--- /dev/null
+++ b/cxxfem/src/CMakeLists.txt
@@ -0,0 +1,12 @@
+
+FILE(GLOB SRCS *.cpp)
+ADD_LIBRARY(fem SHARED ${SRCS})
+TARGET_LINK_LIBRARIES(fem PUBLIC ${GMSH_LIBRARIES})
+TARGET_INCLUDE_DIRECTORIES(fem PUBLIC ${PROJECT_SOURCE_DIR}/src)
+IF(OpenMP_CXX_FOUND)
+ TARGET_LINK_LIBRARIES(fem PUBLIC OpenMP::OpenMP_CXX)
+ENDIF()
+IF(FEM_USE_MKL)
+ TARGET_LINK_LIBRARIES(fem PUBLIC ${MKL_LIBRARIES})
+ TARGET_INCLUDE_DIRECTORIES(fem PUBLIC ${MKL_INCLUDE_DIRS})
+ENDIF()
\ No newline at end of file
diff --git a/cxxfem/src/OpenMP.h b/cxxfem/src/OpenMP.h
new file mode 100644
index 0000000000000000000000000000000000000000..fd3036120c524c3c14bf7b878323898d052d014e
--- /dev/null
+++ b/cxxfem/src/OpenMP.h
@@ -0,0 +1,37 @@
+#ifndef OPENMP_H
+#define OPENMP_H
+
+#include <iostream>
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace OpenMP
+{
+ inline bool enabled()
+ {
+#ifdef _OPENMP
+ return true;
+#else
+ return false;
+#endif
+ }
+
+ inline int get_max_threads()
+ {
+#ifdef _OPENMP
+ return omp_get_max_threads();
+#else
+ return 1;
+#endif
+ }
+ inline void set_num_threads(int n)
+ {
+#ifdef _OPENMP
+ std::cout << "setting num threads to " << n << '\n';
+ return omp_set_num_threads(n);
+#endif
+ }
+};
+
+#endif
diff --git a/cxxfem/src/fem.cpp b/cxxfem/src/fem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9cfd1d4fbff7c25645f3fc536ea59379554f490
--- /dev/null
+++ b/cxxfem/src/fem.cpp
@@ -0,0 +1,14 @@
+#include "fem.h"
+
+namespace fem
+{
+
+FEM_API int debug_level(int lvl)
+{
+ static int level=0;
+ if(lvl>=0)
+ level = lvl;
+ return level;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/fem.h b/cxxfem/src/fem.h
new file mode 100644
index 0000000000000000000000000000000000000000..fd9c2d69004d41f006773e67257af893f5c18f8c
--- /dev/null
+++ b/cxxfem/src/fem.h
@@ -0,0 +1,70 @@
+#ifndef FEM_H
+#define FEM_H
+
+#if defined(WIN32)
+#ifdef fem_EXPORTS
+#define FEM_API __declspec(dllexport)
+#else
+#define FEM_API __declspec(dllimport)
+#endif
+#else
+#define FEM_API
+#endif
+
+#ifdef _MSC_VER
+#pragma warning(disable : 4251) // DLL/templates non exportes
+//#define NOMINMAX // avoids the definition of "min/max" macro in <minwindef.h> (which could shadow std::max)
+#endif //_MSC_VER
+
+#include "config.h"
+#include <vector>
+#include <string>
+#include <iostream>
+
+namespace fem
+{
+class Problem;
+class Entity;
+class ElementSet;
+class Group;
+class Medium;
+class Dirichlet;
+class Neumann;
+class Solver;
+class Partition;
+class Part;
+class DOF;
+class DOFs;
+class DOFMask;
+class Post;
+class Timer;
+class Extractor;
+} // namespace fem
+
+namespace fem
+{
+
+template <typename T>
+std::ostream &
+operator<<(std::ostream &out, std::vector<T> const &v)
+{
+ out << "[ ";
+ for (auto const &i : v)
+ out << i << " ";
+ out << "]";
+ return out;
+}
+
+template <typename T1, typename T2>
+std::ostream &
+operator<<(std::ostream &out, std::pair<T1, T2> const &v)
+{
+ out << "(" << v.first << "," << v.second << ")";
+ return out;
+}
+
+FEM_API int debug_level(int lvl=-1);
+
+} // namespace fem
+
+#endif // FEM_H
\ No newline at end of file
diff --git a/cxxfem/src/femDOF.cpp b/cxxfem/src/femDOF.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..883bfa19ae6d2bd8bcd6be7484196b110268c768
--- /dev/null
+++ b/cxxfem/src/femDOF.cpp
@@ -0,0 +1,49 @@
+#include "femDOF.h"
+
+namespace fem
+{
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, DOF const &obj)
+{
+ out << obj.dofs.name(obj);
+ return out;
+}
+
+DOFs::DOFs()
+{
+}
+
+void DOFs::add(std::string const &name)
+{
+ size_t idx = names.size();
+ names.push_back(name);
+ idx_by_name[name] = idx;
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, DOFMask const &obj)
+{
+ bool first = true;
+ for (size_t i = 0; i < obj.dofs.size(); ++i)
+ {
+ DOF dof = obj.dofs[i];
+ if (obj & dof)
+ {
+ if (!first)
+ out << '|';
+ out << dof;
+ first = false;
+ }
+ }
+ return out;
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, DOFs const &obj)
+{
+ out << "DOF types (size=" << obj.names.size() << "): " << obj.names << '\n';
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femDOF.h b/cxxfem/src/femDOF.h
new file mode 100644
index 0000000000000000000000000000000000000000..2e2ad4fadbfc8b129bb4b3a49ff32688868c74eb
--- /dev/null
+++ b/cxxfem/src/femDOF.h
@@ -0,0 +1,86 @@
+#ifndef FEMDOF_H
+#define FEMDOF_H
+
+#include "fem.h"
+#include <map>
+
+namespace fem
+{
+
+/// a dof type (an index in the dofs table) e.g. "X"=>0, "Y"=>1, etc.
+
+class FEM_API DOF
+{
+ DOFs const &dofs; ///< set of flags
+ size_t idx; ///< readonly index
+public:
+ size_t index() const { return idx; }
+ friend FEM_API std::ostream &operator<<(std::ostream &out, DOF const &obj);
+
+private:
+ friend class DOFs; // DOF must be created by DOFs
+ friend class DOFMask; // used for accessing dofs (for printing)
+ DOF(DOFs const &_dofs, size_t _idx) : dofs(_dofs), idx(_idx) {}
+};
+
+/// stores all dof types
+/// this is the only class which is able to create a new type of dof
+
+class FEM_API DOFs
+{
+ std::vector<std::string> names;
+ std::map<std::string, size_t> idx_by_name;
+
+public:
+ DOFs();
+ ~DOFs() {}
+ size_t size() const { return names.size(); }
+ std::string name(DOF const &dof) const { return names.at(dof.idx); }
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, DOFs const &obj);
+ void add(std::string const &dofname);
+ DOF operator[](std::string const &dofname) const { return DOF(*this, idx_by_name.at(dofname)); }
+ DOF operator[](char const *dofname) const { return DOF(*this, idx_by_name.at(std::string(dofname))); }
+ DOF operator[](size_t idx) const { return DOF(*this, idx); }
+};
+
+/// a mask of dofs
+/// allows the user to select several dof types using |
+
+class FEM_API DOFMask
+{
+ DOFs const &dofs; ///< set
+ size_t v; ///< mask
+public:
+ /// create a mask from a single dof (implicitly if needed)
+ DOFMask(DOF const &dof) : dofs(dof.dofs), v(1 << dof.index()) {}
+ /// combine masks & other dofs
+ DOFMask operator|(DOFMask const &obj) const
+ {
+ if (&dofs != &obj.dofs)
+ throw std::runtime_error("incompatible sets"); // use assert?
+ return DOFMask(dofs, v | obj.v);
+ }
+ /// check a mask versus another one
+ bool operator&(DOFMask const &obj) const
+ {
+ if (&dofs != &obj.dofs)
+ throw std::runtime_error("incompatible sets"); // use assert?
+ return v & obj.v;
+ }
+ /// pretty print
+ friend FEM_API std::ostream &operator<<(std::ostream &out, DOFMask const &obj);
+
+private:
+ /// rebuild a mask from the result of a boolean operation
+ DOFMask(DOFs const &_dofs, size_t _v) : dofs(_dofs), v(_v) {}
+};
+
+inline DOFMask operator|(DOF const &d1, DOF const &d2)
+{
+ return (DOFMask)d1 | d2;
+}
+
+} // namespace fem
+
+#endif // FEMDOF_H
diff --git a/cxxfem/src/femDirichlet.cpp b/cxxfem/src/femDirichlet.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..db33313b6571cabfd004e8e3448b609e209f757d
--- /dev/null
+++ b/cxxfem/src/femDirichlet.cpp
@@ -0,0 +1,29 @@
+#include "femDirichlet.h"
+#include "femProblem.h"
+#include "femGroup.h"
+
+namespace fem
+{
+
+Dirichlet::Dirichlet(Problem &pbl, std::string const &groupname,
+ std::string const &dof, double val) : m_dofname(dof), m_val(val)
+{
+ group = pbl.groups_by_name.at(groupname);
+ pbl.dirichlets.push_back(this);
+}
+
+Dirichlet::~Dirichlet()
+{
+ if(debug_level()>0)
+ std::cout << "~Dirichlet()\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Dirichlet const &obj)
+{
+ out << "Dirichlet BC on \"" << obj.group->name << "\": "
+ << obj.m_dofname << '=' << obj.m_val;
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femDirichlet.h b/cxxfem/src/femDirichlet.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae6e80fd4a3fa37702d529770c2ff75f9c3c5426
--- /dev/null
+++ b/cxxfem/src/femDirichlet.h
@@ -0,0 +1,32 @@
+#ifndef FEMDIRICHLET_H
+#define FEMDIRICHLET_H
+
+#include "fem.h"
+
+namespace fem
+{
+
+/// Dirichlet boundary condition:
+/// a target group and a value to be applied on a given dof
+
+class FEM_API Dirichlet
+{
+public:
+#ifndef SWIG
+ Group *group; ///< target group
+ std::string m_dofname; ///< target degree of freedom
+ double m_val; ///< prescribed value
+#endif
+
+public:
+ Dirichlet(Problem &pbl, std::string const &groupname, std::string const &dof, double val);
+ ~Dirichlet();
+
+#ifndef SWIG
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Dirichlet const &obj);
+#endif
+};
+
+} // namespace fem
+
+#endif // FEMDIRICHLET_H
diff --git a/cxxfem/src/femElementSet.cpp b/cxxfem/src/femElementSet.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd1f0958a8cc01e289cac47c18394db2e4a69c71
--- /dev/null
+++ b/cxxfem/src/femElementSet.cpp
@@ -0,0 +1,48 @@
+#include "femElementSet.h"
+#include "femGmsh.h"
+
+namespace fem
+{
+
+ElementSet::ElementSet(int _typ,
+ std::vector<std::size_t> const &_etags,
+ std::vector<std::size_t> const &_ntags) : type_id(_typ),
+ etags(_etags),
+ ntags_v(_ntags),
+ ntags(&ntags_v[0], etags.size(), ntags_v.size() / etags.size())
+{
+
+ // properties of type "type_id" (could be shared)
+ std::vector<double> localNodeCoord;
+ gmsh::model::mesh::getElementProperties(type_id,
+ type_name, type_dim,
+ type_order, type_nnods,
+ localNodeCoord,
+ type_numPrimaryNodes);
+
+ // localNodeCoord is [xi1, eta1, xi2, eta2, xi3, eta3, ...] for 2d elements
+ // which is inconsistent with localCoords returned for the Gauss points
+ // [xi1, eta1, 0, xi2, eta2, 0, xi3, eta3, 0, ...]
+ type_localNodeCoord.resize(type_nnods * 3);
+ std::fill(type_localNodeCoord.begin(), type_localNodeCoord.end(), 0.0);
+
+ for (int i = 0; i < type_nnods; ++i)
+ for (int j = 0; j < type_dim; ++j)
+ type_localNodeCoord[3 * i + j] = localNodeCoord[type_dim * i + j];
+}
+
+ElementSet::~ElementSet()
+{
+ // std::cout << "~ElementSet()\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, ElementSet const &obj)
+{
+ out << "set of " << obj.etags.size() << " \"" << obj.type_name
+ << "\" (type=" << obj.type_id << ", order=" << obj.type_order
+ << ", nnods=" << obj.type_nnods << ")";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femElementSet.h b/cxxfem/src/femElementSet.h
new file mode 100644
index 0000000000000000000000000000000000000000..878edaa3ac821c090afe71b99d2fcf1d88d93be4
--- /dev/null
+++ b/cxxfem/src/femElementSet.h
@@ -0,0 +1,40 @@
+#ifndef FEMELEMENTSET_H
+#define FEMELEMENTSET_H
+
+#include "fem.h"
+#include <Eigen/Core>
+
+namespace fem
+{
+
+/// a set of elements of a given type
+
+class FEM_API ElementSet
+{
+public:
+ int type_id; ///< type (integer) e.g. 5 => "Hexahedron 8"
+ std::vector<std::size_t> etags; ///< element tags
+ std::vector<std::size_t> ntags_v; ///< node tags [nodes of elem1, then elem2, elem3, etc.]
+
+ Eigen::Map<Eigen::Matrix<size_t, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+ ntags;
+
+ // element type data (could be shared)
+ std::string type_name;
+ int type_dim;
+ int type_order;
+ int type_nnods;
+ std::vector<double> type_localNodeCoord;
+ int type_numPrimaryNodes;
+
+public:
+ ElementSet(int _typ, std::vector<std::size_t> const &_etags,
+ std::vector<std::size_t> const &_ntags);
+ ~ElementSet();
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, ElementSet const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMELEMENTSET_H
diff --git a/cxxfem/src/femEntity.cpp b/cxxfem/src/femEntity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b71cd63adc19c98c70626298f39721fc348e3bb1
--- /dev/null
+++ b/cxxfem/src/femEntity.cpp
@@ -0,0 +1,55 @@
+#include "femEntity.h"
+#include "femElementSet.h"
+#include "femGmsh.h"
+
+namespace fem
+{
+
+Entity::Entity(int _dim, int _tag) : dim(_dim), tag(_tag)
+{
+ gmsh::model::getEntityName(dim, tag, name);
+ gmsh::model::getType(dim, tag, typ);
+
+ // mesh
+ std::vector<int> elementTypes;
+ std::vector<std::vector<std::size_t>> elementTags;
+ std::vector<std::vector<std::size_t>> nodeTags;
+ gmsh::model::mesh::getElements(elementTypes, elementTags, nodeTags, dim, tag);
+ for (size_t i = 0; i < elementTypes.size(); ++i)
+ esets.push_back(new ElementSet(elementTypes[i], elementTags[i], nodeTags[i]));
+}
+
+Entity::~Entity()
+{
+ // std::cout << "~Entity()\n";
+ for (auto &e : esets)
+ delete e;
+}
+
+void Entity::link_to_groups(std::map<std::pair<int, int>, Group *> const &groups_by_dimtag)
+{
+ std::vector<int> tags;
+ gmsh::model::getPhysicalGroupsForEntity(dim, tag, tags);
+ for (auto tag : tags)
+ {
+ try
+ {
+ groups.push_back(groups_by_dimtag.at(std::make_pair(dim, tag)));
+ }
+ catch (const std::exception &e)
+ {
+ throw std::runtime_error("Entity::link_to_groups: entity (" + std::to_string(dim) + "," + std::to_string(tag) + ") not found!");
+ }
+ }
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Entity const &obj)
+{
+ out << obj.typ << " (" << obj.dim << "d, " << obj.tag << ") named \""
+ << obj.name << "\" with " << obj.esets.size() << " mesh(es), belonging to "
+ << obj.groups.size() << " group(s)";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femEntity.h b/cxxfem/src/femEntity.h
new file mode 100644
index 0000000000000000000000000000000000000000..f63cb958be1f62c4d4fe15b4eff2ac9ee0687e7a
--- /dev/null
+++ b/cxxfem/src/femEntity.h
@@ -0,0 +1,33 @@
+#ifndef FEMENTITY_H
+#define FEMENTITY_H
+
+#include "fem.h"
+#include <map>
+
+namespace fem
+{
+
+/// a gmsh entity
+
+class FEM_API Entity
+{
+public:
+ int dim; ///< dimension 0, 1, 2 or 3
+ int tag; ///< gmsh id
+ std::string name; ///< name (usually "")
+ std::string typ; ///< Point, Surface, Plane, Line, Volume, etc.
+ std::vector<ElementSet *> esets; ///< maillages (1 maillage = 1 type d'élément)
+ std::vector<Group *> groups; ///< physical groups
+
+public:
+ Entity(int _dim, int _tag);
+ ~Entity();
+
+ void link_to_groups(std::map<std::pair<int, int>, Group *> const &groups_by_dimtag);
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Entity const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMENTITY_H
diff --git a/cxxfem/src/femExtractor.cpp b/cxxfem/src/femExtractor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fc0a89749bf9ed6aea663fdd2e2d860fe6667bf4
--- /dev/null
+++ b/cxxfem/src/femExtractor.cpp
@@ -0,0 +1,23 @@
+#include "femExtractor.h"
+
+namespace fem
+{
+
+Extractor::Extractor()
+{
+}
+
+Extractor::~Extractor()
+{
+ if(debug_level()>0)
+ std::cout << "~Extractor()\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Extractor const &obj)
+{
+ out << "Extractor\n";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femExtractor.h b/cxxfem/src/femExtractor.h
new file mode 100644
index 0000000000000000000000000000000000000000..6f4348c5cd8544468ec39551fc8c094620210d9f
--- /dev/null
+++ b/cxxfem/src/femExtractor.h
@@ -0,0 +1,23 @@
+#ifndef FEMEXTRACTOR_H
+#define FEMEXTRACTOR_H
+
+#include "fem.h"
+
+namespace fem
+{
+
+///
+
+class FEM_API Extractor
+{
+public:
+public:
+ Extractor();
+ ~Extractor();
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Extractor const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMEXTRACTOR_H
diff --git a/cxxfem/src/femGmsh.h b/cxxfem/src/femGmsh.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6e0428900c06a918547acf90152b2bc2ca2a13b
--- /dev/null
+++ b/cxxfem/src/femGmsh.h
@@ -0,0 +1,5 @@
+#ifdef _MSC_VER
+#include <gmsh.h_cwrap>
+#else
+#include <gmsh.h>
+#endif
diff --git a/cxxfem/src/femGroup.cpp b/cxxfem/src/femGroup.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6ca75075cb6b20cb772e18cb66bc286c4ae37779
--- /dev/null
+++ b/cxxfem/src/femGroup.cpp
@@ -0,0 +1,35 @@
+#include "femGroup.h"
+#include "femGmsh.h"
+
+namespace fem
+{
+
+Group::Group(int _dim, int _tag) : dim(_dim), tag(_tag)
+{
+ gmsh::model::getPhysicalName(dim, tag, name);
+ // this-> entities filled later with 'link_to_entities'
+}
+
+Group::~Group()
+{
+ if(debug_level()>0)
+ std::cout << "~Group()\n";
+}
+
+void Group::link_to_entities(std::map<std::pair<int,int>, Entity*> const &entities_by_dimtag)
+{
+ std::vector<int> tags;
+ gmsh::model::getEntitiesForPhysicalGroup(dim, tag, tags);
+ for(auto tag : tags)
+ entities.push_back( entities_by_dimtag.at(std::make_pair(dim,tag)) );
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Group const &obj)
+{
+ out << "Group \"" << obj.name << "\" (" << obj.dim << "d, " << obj.tag
+ << ") made of " << obj.entities.size() << " entities";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femGroup.h b/cxxfem/src/femGroup.h
new file mode 100644
index 0000000000000000000000000000000000000000..f4c3ac7ffd17b5bc63b55a336d48b8f0326b60a7
--- /dev/null
+++ b/cxxfem/src/femGroup.h
@@ -0,0 +1,31 @@
+#ifndef FEMGROUP_H
+#define FEMGROUP_H
+
+#include "fem.h"
+#include <map>
+
+namespace fem
+{
+
+/// a set of entities (physical group)
+
+class FEM_API Group
+{
+public:
+ int dim; ///< dimension 0, 1, 2 or 3
+ int tag; ///< gmsh id
+ std::string name;
+ std::vector<Entity *> entities;
+
+public:
+ Group(int _dim, int _tag);
+ ~Group();
+
+ void link_to_entities(std::map<std::pair<int, int>, Entity *> const &entities_by_dimtag);
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Group const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMGROUP_H
diff --git a/cxxfem/src/femMaps.h b/cxxfem/src/femMaps.h
new file mode 100644
index 0000000000000000000000000000000000000000..ce25ca3d26bd9a30579e70cc71298eadb266a2c0
--- /dev/null
+++ b/cxxfem/src/femMaps.h
@@ -0,0 +1,114 @@
+#ifndef FEMMAPS_H
+#define FEMMAPS_H
+
+#include "fem.h"
+#include <Eigen/Dense>
+
+namespace fem
+{
+
+/// helper class for accessing gmsh jacobian vector as an Eigen Matrix3d
+/// WARNING: this class does not copy the initial vector!
+///
+/// with Eigen, the jacobian is NOT transposed by default
+/// it is returned "by column" by gmsh and Eigen uses column-major
+/// storage by default
+/// (!= C storage or numpy's default = row major)
+///
+/// gmsh = [ element1, element2, element3, ... ]
+/// element1 = [ gp1, gp2, gp3, ... ]
+/// gp1 = [ dx/du, dy/du, dz/du, dx/dv, ... ]
+///
+/// [ dx/du dx/dv dx/dw ]
+/// JacoMap(e, gp) = [ dy/du dy/dv dy/dw ]
+/// [ dz/du dz/dv dz/dw ]
+///
+class JacoMap
+{
+ std::vector<double> &jacobians;
+ size_t nelems; ///< only used to check bounds in debug - TODO: could be computed from jacobians.size()
+ size_t nGP;
+
+public:
+ typedef Eigen::Map<Eigen::Matrix3d> M;
+
+ JacoMap(std::vector<double> &_jacobians,
+ size_t _nelems, size_t _nGP) : jacobians(_jacobians),
+ nelems(_nelems), nGP(_nGP)
+ {
+ }
+ M operator()(size_t e, size_t gp) const
+ {
+ assert(e < nelems && gp < nGP);
+ return M(&jacobians[(9 * nGP) * e + 9 * gp]);
+ }
+};
+
+/// helper class for accessing gmsh coord vector as an Eigen Vector3d
+/// WARNING: this class does not copy the initial vector!
+///
+/// gmsh = [ elem1, elem2, elem3, ... ]
+/// elem1 = [ gp1, gp2, gp3, ... ]
+/// gp1 = [ x, y, z ]
+///
+/// CoordMap(e, gp) = [ x ]
+/// [ y ]
+/// [ z ]
+
+class CoordMap
+{
+ std::vector<double> &coord;
+ size_t nelems; ///< only used to check bounds in debug - TODO: could be computed from coord.size()
+ size_t nGP;
+
+public:
+ typedef Eigen::Map<Eigen::Vector3d> M;
+
+ CoordMap(std::vector<double> &_coord,
+ size_t _nelems, size_t _nGP) : coord(_coord),
+ nelems(_nelems), nGP(_nGP)
+ {
+ }
+ M operator()(size_t e, size_t gp) const
+ {
+ assert(e < nelems && gp < nGP);
+ return M(&coord[(3 * nGP) * e + 3 * gp]);
+ }
+};
+
+/// helper class for accessing gmsh dshapef vector as an Eigen MatrixXd
+///
+/// gmsh = [ gp1, gp2, gp2, ... ]
+/// gp1 = [ shp1, shp2, shp3, ... ]
+/// shp1 = [df/du, df/dv, df/dw]
+///
+/// DShapeMap(gp) = [ df/du, df/dv, df/dw ] <= shp1
+/// [ df/du, df/dv, df/dw ] <= shp2
+/// [ df/du, df/dv, df/dw ] <= shp3
+/// [ ... ]
+
+class DShapeMap
+{
+ std::vector<double> &basisDF;
+ size_t nGP; ///< only used to check bounds in debug - TODO: could be computed from basisDF.size()
+ int nnod;
+
+public:
+ typedef Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> M;
+
+ DShapeMap(std::vector<double> &_basisDF,
+ size_t _nGP, int _nnod) : basisDF(_basisDF), nGP(_nGP),
+ nnod(_nnod)
+ {
+ }
+
+ M operator()(size_t gp) const
+ {
+ assert(gp < nGP);
+ return M(&basisDF[(nnod * 3) * gp], nnod, 3);
+ }
+};
+
+} // namespace fem
+
+#endif // FEMMAPS_H
diff --git a/cxxfem/src/femMedium.cpp b/cxxfem/src/femMedium.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d6037909cbecc6bd792b33cf3eaccc79c02c086a
--- /dev/null
+++ b/cxxfem/src/femMedium.cpp
@@ -0,0 +1,50 @@
+#include "femMedium.h"
+#include "femProblem.h"
+#include "femGroup.h"
+
+namespace fem
+{
+
+Medium::Medium(Problem &pbl, std::string const &groupname,
+ double E, double nu) : m_E(E), m_nu(nu), H2D(3, 3), H3D(6, 6)
+{
+ if (pbl.entities.empty()) // is the problem already sync'ed?
+ pbl.sync();
+
+ group = pbl.groups_by_name.at(groupname);
+ pbl.media.push_back(this);
+
+ // construct 2D-plane-stress Hooke matrix
+ H2D << 1.0, nu, 0.0,
+ nu, 1, 0.0,
+ 0.0, 0.0, (1.0 - nu) / 2.0;
+ H2D *= E / (1 - nu * nu);
+
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "H2D =" << H2D.format(fmt);
+
+ // construct 3D Hooke matrix
+ H3D << 1.0 - nu, nu, nu, 0.0, 0.0, 0.0,
+ nu, 1.0 - nu, nu, 0.0, 0.0, 0.0,
+ nu, nu, 1.0 - nu, 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, (1.0 - 2 * nu) / 2.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, (1.0 - 2 * nu) / 2.0, 0.0,
+ 0.0, 0.0, 0.0, 0.0, 0.0, (1.0 - 2 * nu) / 2.0;
+ H3D *= E / (1 + nu) / (1 - 2 * nu);
+}
+
+Medium::~Medium()
+{
+ if(debug_level()>0)
+ std::cout << "~Medium()\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Medium const &obj)
+{
+ out << "Medium on \"" << obj.group->name
+ << "\": E=" << obj.m_E << ", nu=" << obj.m_nu;
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femMedium.h b/cxxfem/src/femMedium.h
new file mode 100644
index 0000000000000000000000000000000000000000..88e779cce05787dec1d06b6b6cafb2f1b546b186
--- /dev/null
+++ b/cxxfem/src/femMedium.h
@@ -0,0 +1,40 @@
+#ifndef FEMMEDIUM_H
+#define FEMMEDIUM_H
+
+#include "fem.h"
+#include <Eigen/Dense>
+
+namespace fem
+{
+
+/// physical properties on a domain defined by a group
+
+class FEM_API Medium
+{
+public:
+#ifndef SWIG
+ double m_E; ///< Young's modulus
+ double m_nu; ///< Poisson's ratio
+ Group *group;
+
+ Eigen::MatrixXd H2D; ///< Hooke
+ Eigen::MatrixXd H3D; ///< Hooke
+#endif
+
+public:
+ Medium(Problem &pbl, std::string const &groupname, double E, double nu);
+ ~Medium();
+
+#ifndef SWIG
+ Eigen::MatrixXd &hooke(int ndim)
+ {
+ return (ndim==2)? H2D : H3D;
+ }
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Medium const &obj);
+#endif
+};
+
+} // namespace fem
+
+#endif // FEMMEDIUM_H
diff --git a/cxxfem/src/femNeumann.cpp b/cxxfem/src/femNeumann.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2a2b6948fdeefb2be4b1abbabbadfd0d6fc1c1e9
--- /dev/null
+++ b/cxxfem/src/femNeumann.cpp
@@ -0,0 +1,49 @@
+#include "femNeumann.h"
+#include "femProblem.h"
+#include "femGroup.h"
+
+namespace fem
+{
+
+Neumann::Neumann(Problem &pbl, std::string const &groupname,
+ std::string const &dof, double val)
+ : m_nodetag(0), m_dofname(dof), m_val(val)
+{
+ try
+ {
+ group = pbl.groups_by_name.at(groupname);
+ }
+ catch (const std::exception &e)
+ {
+ throw std::runtime_error("Neumann: group \"" + groupname + "\" not found!");
+ }
+ pbl.neumanns.push_back(this);
+}
+
+Neumann::Neumann(Problem &pbl, size_t nodetag,
+ std::string const &dof, double val)
+ : group(nullptr), m_nodetag(nodetag), m_dofname(dof), m_val(val)
+{
+ pbl.neumanns.push_back(this);
+}
+
+Neumann::~Neumann()
+{
+ // if(debug_level()>0)
+ // std::cout << "~Neumann()\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Neumann const &obj)
+{
+ if (obj.group)
+ out << "Neumann BC on \"" << obj.group->name << "\": "
+ << obj.m_dofname << '=' << obj.m_val;
+ // else
+ // out << "Neumann BC on node " << obj.m_nodetag << ": "
+ // << obj.m_dofname << '=' << obj.m_val;
+
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femNeumann.h b/cxxfem/src/femNeumann.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a51cad5900a6a3f1e78961b9795e49f9f468603
--- /dev/null
+++ b/cxxfem/src/femNeumann.h
@@ -0,0 +1,34 @@
+#ifndef FEMNEUMANN_H
+#define FEMNEUMANN_H
+
+#include "fem.h"
+
+namespace fem
+{
+
+/// Neumann boundary condition:
+/// a target group and a value to be applied on a given dof
+
+class FEM_API Neumann
+{
+public:
+#ifndef SWIG
+ Group *group; ///< option1: target group
+ size_t m_nodetag; ///< option2: target node tag
+ std::string m_dofname; ///< target degree of freedom
+ double m_val; ///< prescribed value
+#endif
+
+public:
+ Neumann(Problem &pbl, std::string const &groupname, std::string const &dof, double val);
+ Neumann(Problem &pbl, size_t nodetag, std::string const &dof, double val);
+ ~Neumann();
+
+#ifndef SWIG
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Neumann const &obj);
+#endif
+};
+
+} // namespace fem
+
+#endif // FEMNEUMANN_H
diff --git a/cxxfem/src/femNewClass.cpp b/cxxfem/src/femNewClass.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..22db923ca26a0a5100492428a025437f73225e63
--- /dev/null
+++ b/cxxfem/src/femNewClass.cpp
@@ -0,0 +1,23 @@
+#include "femNewClass.h"
+
+namespace fem
+{
+
+NewClass::NewClass()
+{
+}
+
+NewClass::~NewClass()
+{
+ if(debug_level()>0)
+ std::cout << "~NewClass()\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, NewClass const &obj)
+{
+ out << "NewClass\n";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femNewClass.h b/cxxfem/src/femNewClass.h
new file mode 100644
index 0000000000000000000000000000000000000000..4decb13d9db89f6a333474fa49b6f0978476e64c
--- /dev/null
+++ b/cxxfem/src/femNewClass.h
@@ -0,0 +1,23 @@
+#ifndef FEMNEWCLASS_H
+#define FEMNEWCLASS_H
+
+#include "fem.h"
+
+namespace fem
+{
+
+///
+
+class FEM_API NewClass
+{
+public:
+public:
+ NewClass();
+ ~NewClass();
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, NewClass const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMNEWCLASS_H
diff --git a/cxxfem/src/femPart.cpp b/cxxfem/src/femPart.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b872b8d21e997a14f3654228e82fdfae7a151118
--- /dev/null
+++ b/cxxfem/src/femPart.cpp
@@ -0,0 +1,28 @@
+#include "femPart.h"
+#include "femPartition.h"
+
+namespace fem
+{
+
+Part::Part(Partition &partition,
+ DOFMask const &_dofmask,
+ DOFMask const &_bcmask) : dofmask(_dofmask), bcmask(_bcmask)
+{
+ idx = partition.parts.size();
+ partition.parts.push_back(this);
+}
+
+Part::~Part()
+{
+ if(debug_level()>0)
+ std::cout << "~Part()\n";
+}
+
+std::ostream &
+operator<<(std::ostream &out, Part const &obj)
+{
+ out << "Part #" << obj.idx << ": " << obj.dofmask << ", " << obj.bcmask;
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femPart.h b/cxxfem/src/femPart.h
new file mode 100644
index 0000000000000000000000000000000000000000..79555e147e49b14d929b6318cc1452b204c3a5ea
--- /dev/null
+++ b/cxxfem/src/femPart.h
@@ -0,0 +1,28 @@
+#ifndef FEMPART_H
+#define FEMPART_H
+
+#include "fem.h"
+#include "femDOF.h"
+
+namespace fem
+{
+
+///
+
+class FEM_API Part
+{
+public:
+ DOFMask dofmask; ///< DOF type
+ DOFMask bcmask; ///< BC type
+ size_t idx; ///< partition number (index in partition.parts)
+
+public:
+ Part(Partition &partition, DOFMask const &_dofmask, DOFMask const &_bcmask);
+ ~Part();
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Part const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMPART_H
diff --git a/cxxfem/src/femPartition.cpp b/cxxfem/src/femPartition.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..90298e9789dc782cb53a4687fe6fbfd407c04eb9
--- /dev/null
+++ b/cxxfem/src/femPartition.cpp
@@ -0,0 +1,26 @@
+#include "femPartition.h"
+#include "femPart.h"
+
+namespace fem
+{
+
+Partition::Partition()
+{
+}
+
+Partition::~Partition()
+{
+ if(debug_level()>0)
+ std::cout << "~Partition()\n";
+}
+
+std::ostream &
+operator<<(std::ostream &out, Partition const &obj)
+{
+ out << "Partition\n";
+ for(auto p : obj.parts)
+ out << '\t' << *p << '\n';
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femPartition.h b/cxxfem/src/femPartition.h
new file mode 100644
index 0000000000000000000000000000000000000000..24c9360f6d0679c5686979ef98beec7032af2172
--- /dev/null
+++ b/cxxfem/src/femPartition.h
@@ -0,0 +1,25 @@
+#ifndef FEMPARTITION_H
+#define FEMPARTITION_H
+
+#include "fem.h"
+
+namespace fem
+{
+
+/// a partition is a set of parts
+
+class FEM_API Partition
+{
+public:
+ std::vector<Part *> parts;
+
+public:
+ Partition();
+ ~Partition();
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Partition const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMPARTITION_H
diff --git a/cxxfem/src/femPost.cpp b/cxxfem/src/femPost.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a421fb3997faf9f372db711111767a7718179e48
--- /dev/null
+++ b/cxxfem/src/femPost.cpp
@@ -0,0 +1,409 @@
+#include "femPost.h"
+#include "femSolver.h"
+#include "femProblem.h"
+#include "femGmsh.h"
+#include "femGroup.h"
+#include "femEntity.h"
+#include <numeric>
+
+namespace fem
+{
+
+Post::Post(Solver &_solver) : verbose(false), solver(_solver)
+{
+ // global gmsh options
+ gmsh::option::setNumber("General.Orthographic", 0); // better perspective
+ gmsh::option::setNumber("Mesh.VolumeEdges", 0); // don't display the mesh
+ gmsh::option::setNumber("Mesh.SurfaceEdges", 0); // don't display the mesh
+}
+
+Post::~Post()
+{
+ if (debug_level() > 0)
+ std::cout << "~Post()\n";
+}
+
+/// rebuild solution vectors for gmsh and add views
+
+void Post::build_views()
+{
+ std::cout << "creating views...\n";
+ // nodal dofs
+ add_dofs_view(solver.pbl, solver.dofs, solver.X, solver.dofPart, solver.dofIdx);
+
+ // nodal vectors
+ add_nodalvector_view(solver.pbl, "displacement", solver.dofs, solver.X, solver.dofPart, solver.dofIdx);
+ add_nodalvector_view(solver.pbl, "force", solver.dofs, solver.f, solver.dofPart, solver.dofIdx);
+
+ // compute heat fluxes, strains, stresses
+
+ std::vector<size_t> etags_all; // tags of elements
+ std::vector<std::vector<double>> sigma_all; // values of stresses at the nodes of each element
+ std::vector<std::vector<double>> epsilon_all; // values of strains at the nodes of each element
+ std::map<size_t, std::vector<Eigen::VectorXd>> sig_by_node; // value of all the stress vectors (Voigt) at a given node (key of the dict)
+ std::map<size_t, std::vector<Eigen::VectorXd>> eps_by_node; // same for strain vector
+ solver.compute_sig_epl(etags_all, sigma_all, epsilon_all, sig_by_node, eps_by_node);
+
+ // ...discontinuous by node by element
+ add_tensor_views("stress", etags_all, sigma_all);
+ add_tensor_views("strain", etags_all, epsilon_all);
+
+ // ...continuous averaged nodal fields
+ add_smooth_views("stress", sig_by_node);
+ add_smooth_views("strain", eps_by_node);
+}
+
+/// add 1 scalar view for each dof type to gmsh.
+
+void Post::add_dofs_view(Problem &pbl, DOFs &dofs,
+ std::vector<Eigen::VectorXd *> &X,
+ Eigen::ArrayXXi &dofPart,
+ Eigen::ArrayXXi &dofIdx)
+{
+ // scalar view of each dof
+ for (size_t di = 0; di < dofs.size(); ++di) // TODO: iterator
+ {
+ DOF dof = dofs[di]; // extract each dof type separately
+ std::vector<double> x_all(pbl.nodeTags.size(), 0.0);
+ for (size_t i = 0; i < pbl.nodeTags.size(); ++i)
+ {
+ size_t ntag = pbl.nodeTags[i];
+ size_t gidx = pbl.nodemap[ntag];
+ int part = dofPart(gidx, dof.index());
+ int idx = dofIdx(gidx, dof.index());
+ x_all[i] = (*X[part])[idx];
+ }
+ // Export solution to gmsh
+ std::string name = dofs.name(dof);
+ int sview = add_view(name);
+ gmsh::view::addHomogeneousModelData(sview, 0, "", "NodeData",
+ pbl.nodeTags, x_all);
+ }
+}
+
+// add 1 tensor view to gmsh and 6 scalar views for the components.
+// tensor_all is in gmsh format: array of "nbelem" elements (same size as "etags_all")
+// each element contains the tensor data:
+// [a11 a12 a13 a21 a22 a23 a31 a32 a33 a11 a12 a13 a21 a22 a23 a31 a32 a33 ...]
+// <--------------node1--------------> <--------------node2--------------> <--
+
+void Post::add_tensor_views(std::string const &vname,
+ std::vector<size_t> const &etags_all,
+ std::vector<std::vector<double>> const &tensor_all)
+{
+ // std::cout << "Post::add_tensor_views " << vname << "\n";
+ // std::cout << "\tetags_all.size() = " << etags_all.size() << "\n";
+ // std::cout << "\ttensor_all.size() = " << tensor_all.size() << "\n";
+
+ // tensor view (by node, by element)
+ std::string name = vname + "_tensor";
+ int sview = add_view(name);
+ gmsh::view::addModelData(sview, 0, "", "ElementNodeData",
+ etags_all, tensor_all, 0.0, 9);
+
+ // component view (by node, by element)
+ // components are retrieved from the tensor data in gmsh format
+ // [ [0] [ [XX] [ [0] [1] [2] [ [XX] [XY] [XZ]
+ // [4] [YY] [3] [4] [5] [XY] [YY] [YZ]
+ // [8] [ZZ] [6] [7] [8] ] [XZ] [YZ] [ZZ] ]
+ // [1or3] [XY]
+ // [2or6] [XZ]
+ // [5or7] ] [YZ] ]
+
+ std::vector<int> comp_id = {0, 4, 8, 1, 2, 5};
+ std::vector<std::string> comp_name = {"xx", "yy", "zz", "xy", "xz", "yz"};
+ for (size_t i = 0; i < comp_id.size(); ++i)
+ {
+ std::string name = vname + "_" + comp_name[i];
+ int sview = add_view(name);
+ std::vector<std::vector<double>> comp(tensor_all.size());
+ for (size_t ie = 0; ie < tensor_all.size(); ++ie)
+ {
+ comp[ie].resize(tensor_all[ie].size() / 9);
+ for (size_t in = 0; in < tensor_all[ie].size() / 9; ++in)
+ comp[ie][in] = tensor_all[ie][9 * in + comp_id[i]];
+ }
+ gmsh::view::addModelData(sview, 0, "", "ElementNodeData",
+ etags_all, comp, 0.0, 1);
+ }
+}
+
+/// add the nodal vector (displacement, forces) as a new gmsh view.
+
+void Post::add_nodalvector_view(Problem &pbl,
+ std::string const &vname,
+ DOFs &dofs,
+ std::vector<Eigen::VectorXd *> &X,
+ Eigen::ArrayXXi &dofPart,
+ Eigen::ArrayXXi &dofIdx)
+{
+ std::vector<std::vector<double>> uvw_all(pbl.nodeTags.size());
+ for (size_t i = 0; i < pbl.nodeTags.size(); ++i)
+ {
+ uvw_all[i].resize(3);
+ size_t ntag = pbl.nodeTags[i];
+ size_t gidx = pbl.nodemap[ntag];
+ for (size_t di = 0; di < dofs.size(); ++di)
+ {
+ DOF dof = dofs[di];
+ int part = dofPart(gidx, dof.index());
+ int idx = dofIdx(gidx, dof.index());
+ uvw_all[i][di] = (*X[part])[idx];
+ }
+ }
+ // Export solution to gmsh
+ std::string name = vname + "_vector";
+ int sview = add_view(name);
+ gmsh::view::addModelData(sview, 0, "", "NodeData",
+ pbl.nodeTags, uvw_all);
+}
+
+/// add 6 scalar views for the components of a smoothed field (stresses or strains)
+/// voigt_by_node is a map such that:
+/// voigt_by_node[nodetag] = array of vectors of stresses/strains (Voigt - 6 comps)
+/// evaluated from the neighbouring elements
+
+void Post::add_smooth_views(std::string const &vname,
+ std::map<size_t, std::vector<Eigen::VectorXd>> &voigt_by_node)
+{
+ std::vector<Eigen::VectorXd> smooth_sig;
+ std::vector<size_t> smooth_ntags;
+
+ std::vector<std::vector<double>> tensors;
+
+ // compute the average
+
+ for (auto const &[ntag, sigs] : voigt_by_node)
+ {
+ // auto s = std::reduce(sigs.begin(), sigs.end()) / sigs.size(); // C++17 / not avail in gcc 8.1 mingw
+ Eigen::VectorXd s0 = Eigen::VectorXd::Zero(6);
+ // std::cout << "ntag=" << ntag << " sigs.size() = " << sigs.size() << std::endl;
+ // auto s = std::accumulate(sigs.begin(), sigs.end(), s0) / sigs.size(); // core dump sous gcc 9.4 ubuntu ???
+ for (auto &v : sigs)
+ s0 += v;
+ auto s = s0 / sigs.size();
+ // std::cout << "s=" << s << std::endl;
+ smooth_sig.push_back(s);
+ smooth_ntags.push_back(ntag);
+
+ tensors.push_back({s[0], s[3], s[4], s[3], s[1], s[5], s[4], s[5], s[2]});
+ }
+
+ // tensor view (by node)
+ std::string name = "smooth_" + vname + "_tensor";
+ int sview = add_view(name);
+ gmsh::view::addModelData(sview, 0, "", "NodeData",
+ smooth_ntags, tensors, 0.0, 9);
+
+ // loop over the 6 components
+
+ std::vector<std::string> comps = {"xx", "yy", "zz", "xy", "xz", "yz"};
+
+ for (size_t c = 0; c < comps.size(); ++c)
+ {
+ std::string name = "smooth_" + vname + "_" + comps[c];
+ int sview = add_view(name);
+ std::vector<double> comp;
+ for (auto &s : smooth_sig)
+ comp.push_back(s[c]);
+ gmsh::view::addHomogeneousModelData(sview, 0, "", "NodeData",
+ smooth_ntags, comp);
+ }
+
+ // TODO: we could rebuild the average tensor here [DONE]
+ // we could also calculate the 2nd invariant (Von Mises)
+}
+
+/// add a new gmsh view and configure it
+/// store the view in a dict called 'views' so that the views can be easily configured
+/// later by the user using their names
+
+int Post::add_view(std::string const &name)
+{
+ // if(self.verbose): print(f'\tadding new view: {name}')
+
+ int sview = gmsh::view::add(name);
+ // sview = 1 for the first view although it is named Views[0] in the GUI and the options!!
+ views[name] = sview;
+ configure(sview);
+ return sview;
+}
+
+void Post::set_camera_3D()
+{
+ gmsh::option::setNumber("General.TrackballQuaternion0", 0.62);
+ gmsh::option::setNumber("General.TrackballQuaternion1", 0.17);
+ gmsh::option::setNumber("General.TrackballQuaternion2", 0.20);
+ gmsh::option::setNumber("General.TrackballQuaternion3", 0.74);
+}
+
+void Post::deform(double factor)
+{
+ int v_disp = views["displacement_vector"];
+ for (auto const &[vname, v] : views)
+ {
+ if (v != v_disp)
+ {
+ gmsh::option::setNumber("View[" + std::to_string(v - 1) + "].UseGeneralizedRaise", 1);
+ gmsh::option::setNumber("View[" + std::to_string(v - 1) + "].GeneralizedRaiseView", v_disp - 1);
+ gmsh::option::setNumber("View[" + std::to_string(v - 1) + "].GeneralizedRaiseFactor", factor);
+ }
+ }
+
+ // v_disp = self.views['displacement_vector']
+ // for vname in self.views:
+ // v = self.views[vname]
+ // if v != v_disp:
+ // gmsh.option.setNumber(f"View[{v-1}].UseGeneralizedRaise", 1)
+ // gmsh.option.setNumber(f"View[{v-1}].GeneralizedRaiseView", v_disp-1)
+ // gmsh.option.setNumber(f"View[{v-1}].GeneralizedRaiseFactor", factor)
+}
+
+void Post::configure(int sview)
+{
+ // scalar views
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].IntervalsType", 3); // filled iso
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].NbIso", 20); //
+
+ // vector views
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].GlyphLocation", 2); // at the node (default=barycenter)
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].VectorType", 2); // simple Arrow (default=3D arrow)
+
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].ShowElement", 1); // view black mesh outline
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].DrawSkinOnly", 1); // display skin only (should be faster)
+ gmsh::option::setNumber("View[" + std::to_string(sview - 1) + "].Visible", 0); // hide all views
+}
+
+/// display results
+
+void Post::view()
+{
+ std::cout << "viewing results...\n";
+ gmsh::fltk::run();
+}
+
+void Post::show_views(std::vector<std::string> const &names)
+{
+ for (auto &v : names)
+ {
+ auto it = views.find(v);
+ if (it != views.end())
+ gmsh::option::setNumber("View[" + std::to_string(it->second - 1) + "].Visible", 1);
+ else
+ std::cout << "WARNING: view named \"" << v << "\" does not exist!\n";
+ }
+}
+
+/// write mesh/data to disk
+
+void Post::write()
+{
+ gmsh::option::setNumber("Mesh.Binary", 1); // writes files in binary format
+
+ gmsh::write("mesh.opt"); // save the options
+ gmsh::write("mesh.msh"); // only save the mesh
+ for (auto const &[name, sview] : views)
+ {
+ // gmsh::view::write(sview, name+".pos"); // pos format: list-based data (indep. of any model)
+ gmsh::view::write(sview, name + ".msh"); // writes the mesh in each file - each view is in a separate file
+ // gmsh::view::write(sview, "mesh.msh", true); // single file, but the mesh is repeated
+ }
+}
+
+std::vector<double>
+Post::probe(std::string const &name, std::string const &grpname)
+{
+ // std::cout << "Post::probe(" << grpname << ")\n";
+ // get physical group from grpname
+ Group *grp = solver.pbl.groups_by_name.at(grpname); // TODO: better handling of errors
+
+ // get the view
+ auto it = views.find(name);
+ if (it == views.end())
+ throw std::runtime_error("WARNING: view named \"" + name + "\" does not exist!");
+ int nview = it->second;
+
+ // probe values
+ std::vector<double> all_values;
+
+ std::vector<std::size_t> nodeTags;
+ std::vector<double> coord;
+ gmsh::model::mesh::getNodesForPhysicalGroup(grp->dim, grp->tag, nodeTags, coord);
+ // j'ai essayé de boucler sur les entités du groupe et faire un getNode de l'entité
+ // mais ça ne me retourne rien (?)
+
+ // std::cout << "nodeTags.size()=" << nodeTags.size() << '\n';
+
+ for (size_t i = 0; i < nodeTags.size(); ++i)
+ {
+ std::vector<double> values;
+ double distance = 0.0;
+ int step = -1;
+ int numComp = -1; // do not use numComp=comp (gmsh returns [ ] for comp!=0 )
+ bool gradient = false;
+ double distanceMax = -1.; // find closest node
+ gmsh::view::probe(nview, coord[i * 3], coord[i * 3 + 1], coord[i * 3 + 2], values,
+ distance, step, numComp, gradient, distanceMax);
+
+ // debug
+ // std::cout << "probing " << name << " at node #" << nodeTags[i] << ":";
+ // std::cout << "\tnview: " << nview << '\n';
+ // std::cout << "\tcoords: " << coord[i * 3] << ", " << coord[i * 3 + 1] << ", " << coord[i * 3 + 2] << '\n';
+ // std::cout << "\tgot from gmsh: " << values << '\n';
+ // std::cout << "\tdistance = " << distance << std::endl;
+
+ all_values.insert(all_values.end(), values.begin(), values.end());
+ }
+
+ return all_values;
+}
+
+/// Probe a field at a given nodetag
+
+std::vector<double>
+Post::probe(std::string const &name, size_t nodetag)
+{
+ // get the view
+ auto it = views.find(name);
+ if (it == views.end())
+ throw std::runtime_error("WARNING: view named \"" + name + "\" does not exist!");
+ int nview = it->second;
+
+ // retrieve node coordinates
+ std::vector<double> coord;
+ std::vector<double> parametricCoord; // unused
+ int dim; // unused
+ int tag; // unused
+ gmsh::model::mesh::getNode(nodetag, coord, parametricCoord, dim, tag);
+
+ // probe value
+
+ std::vector<double> values;
+ double distance = 0.0;
+ int step = -1;
+ int numComp = -1; // do not use numComp=comp (gmsh returns [ ] for comp!=0 )
+ bool gradient = false;
+ double distanceMax = -1.; // find closest node
+ gmsh::view::probe(nview, coord[0], coord[1], coord[2], values, distance, step, numComp, gradient, distanceMax);
+
+ // debug
+ // std::cout << "probing " << name << " at node #" << nodetag << ":";
+ // std::cout << "\tnview: "<< nview << '\n';
+ // std::cout << "\tcoords: "<< coord << '\n';
+ // std::cout << "\tgot from gmsh: "<< values << '\n';
+ // std::cout << "\tdistance = "<< distance << std::endl;
+
+ // std::cout << "probing \"" << name << "\" (comp=" << comp
+ // << ") at node #" << nodetag << ": " << values[comp] << std::endl;
+
+ return values;
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Post const &obj)
+{
+ out << "Post\n";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femPost.h b/cxxfem/src/femPost.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f48d93f831f1fc8731c9078a85c6342762ed243
--- /dev/null
+++ b/cxxfem/src/femPost.h
@@ -0,0 +1,64 @@
+#ifndef FEMPOST_H
+#define FEMPOST_H
+
+#include "fem.h"
+#include <map>
+#include <list>
+#include <Eigen/Dense>
+
+namespace fem
+{
+
+///
+
+class FEM_API Post
+{
+public:
+#ifndef SWIG
+ bool verbose;
+ Solver &solver;
+ std::map<std::string, int> views;
+#endif
+
+public:
+ Post(Solver &_solver);
+ ~Post();
+
+ void build_views();
+ void view();
+ void show_views(std::vector<std::string> const &names);
+ void set_camera_3D();
+ void deform(double factor=1.0);
+ void write();
+
+ std::vector<double> probe(std::string const &name, std::string const &grpname);
+ std::vector<double> probe(std::string const &name, size_t nodetag);
+
+#ifndef SWIG
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Post const &obj);
+#endif
+
+private:
+ void add_dofs_view(Problem &pbl, DOFs &dofs,
+ std::vector<Eigen::VectorXd *> &X,
+ Eigen::ArrayXXi &dofPart,
+ Eigen::ArrayXXi &dofIdx);
+ void add_nodalvector_view(Problem &pbl,
+ std::string const &vname,
+ DOFs &dofs,
+ std::vector<Eigen::VectorXd *> &X,
+ Eigen::ArrayXXi &dofPart,
+ Eigen::ArrayXXi &dofIdx);
+ void add_tensor_views(std::string const &vname,
+ std::vector<size_t> const &etags_all,
+ std::vector<std::vector<double>> const &tensor_all);
+ void add_smooth_views(std::string const &vname,
+ std::map<size_t, std::vector<Eigen::VectorXd>> &voigt_by_node);
+
+ int add_view(std::string const &name);
+ void configure(int sview);
+};
+
+} // namespace fem
+
+#endif // FEMPOST_H
diff --git a/cxxfem/src/femProblem.cpp b/cxxfem/src/femProblem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d120a18cf008cc58a34436b8583b1ee724f5695
--- /dev/null
+++ b/cxxfem/src/femProblem.cpp
@@ -0,0 +1,148 @@
+#include "femProblem.h"
+#include "femEntity.h"
+#include "femElementSet.h"
+#include "femGroup.h"
+#include "femMedium.h"
+#include "femDirichlet.h"
+#include "femNeumann.h"
+#include "femGmsh.h"
+#include <fstream>
+
+namespace fem
+{
+
+Problem::Problem() : quadrature("Gauss2")
+{
+ gmsh::initialize();
+ // General.Verbosity
+ // (0: silent except for fatal errors, 1: +errors, 2: +warnings, 3: +direct, 4: +information, 5: +status (default), 99: +debug)
+ gmsh::option::setNumber("General.Verbosity", 2);
+ gmsh::option::setNumber("General.Terminal", 1); // Should information be printed on the terminal? 0: no (default) ; 1: yes
+}
+
+Problem::~Problem()
+{
+ for (auto &e : entities)
+ delete e;
+}
+
+/// fill internal helper objects from gmsh data
+
+void Problem::sync()
+{
+ // -- entities
+ std::cout << "retrieving entities...\n";
+ gmsh::vectorpair dimTags;
+ gmsh::model::getEntities(dimTags);
+ for (auto &p : dimTags)
+ entities.push_back(new Entity(p.first, p.second));
+ // create useful maps
+ for (auto e : entities)
+ entities_by_dimtag[std::make_pair(e->dim, e->tag)] = e;
+
+ // -- physical groups
+ std::cout << "retrieving physical groups...\n";
+ dimTags.clear();
+ gmsh::model::getPhysicalGroups(dimTags);
+ for (auto &p : dimTags)
+ groups.push_back(new Group(p.first, p.second));
+ // create useful maps
+ for (auto g : groups)
+ {
+ groups_by_dimtag[std::make_pair(g->dim, g->tag)] = g;
+ groups_by_name[g->name] = g;
+ }
+
+ // link entities <=> groups
+ for (auto g : groups)
+ g->link_to_entities(entities_by_dimtag);
+ for (auto e : entities)
+ e->link_to_groups(groups_by_dimtag);
+
+ // -- nodes
+ // get all the node tags and create a map ("node tag" => "index")
+ // that could be used to store things at nodes in separate arrays
+ std::vector<double> coord;
+ std::vector<double> parametricCoord;
+ gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord);
+ for (size_t i = 0; i < nodeTags.size(); ++i)
+ nodemap[nodeTags[i]] = i;
+
+ //std::cout << *this;
+}
+
+void Problem::load(std::string const &filename,
+ std::map<std::string, double> const &pars,
+ int order)
+{
+ // std::cout << "pars=\n";
+ // for (auto &p : pars)
+ // std::cout << p.first << " : " << p.second << '\n';
+
+ std::string fext;
+ auto idx = filename.rfind('.');
+ if(idx != std::string::npos)
+ fext = filename.substr(idx);
+
+ if (fext!=".geo" && fext!=".msh")
+ throw std::runtime_error("unknown file extension ("+fext+")");
+
+ // load file and initialise gmsh
+ // gmsh::initialize();
+ for (auto &p : pars)
+ {
+ if(debug_level()>1)
+ std::cout << "setting parameter " << p.first << "=" << p.second << '\n';
+ std::vector<double> v;
+ v.push_back(p.second);
+ gmsh::onelab::setNumber(p.first, v);
+ }
+ std::cout << "loading " << filename << '\n';
+ if(! std::ifstream(filename.c_str()).good() ) // not tested by gmsh
+ throw std::runtime_error("file " + filename + " does not exist");
+ gmsh::open(filename);
+
+ // generate mesh
+ if (fext==".geo")
+ {
+ gmsh::model::mesh::generate(3);
+ gmsh::model::mesh::setOrder(order); // after generation
+ quadrature = "Gauss" + std::to_string(2*order); // change default quadrature rule
+ }
+
+ // do not sync now, so that we can add things such as
+ // additional physical groups
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Problem const &obj)
+{
+ out << "Problem:\n";
+ out << '\t' << obj.entities.size() << " entities\n";
+ if(debug_level()>1)
+ for (auto e : obj.entities)
+ {
+ out << "\t\t" << *e << '\n';
+ for (auto m : e->esets)
+ out << "\t\t\tmesh: " << *m << '\n';
+ }
+ out << '\t' << obj.groups.size() << " groups\n";
+ if(debug_level()>1)
+ for (auto g : obj.groups)
+ out << "\t\t" << *g << '\n';
+ out << '\t' << obj.media.size() << " media\n";
+ if(debug_level()>1)
+ for (auto m : obj.media)
+ out << "\t\t" << *m << '\n';
+ out << '\t' << obj.dirichlets.size() << " dirichlets\n";
+ if(debug_level()>1)
+ for (auto bc : obj.dirichlets)
+ out << "\t\t" << *bc << '\n';
+ out << '\t' << obj.neumanns.size() << " neumanns\n";
+ if(debug_level()>1) // TODO: si >20 ... afficher les 10 premiers et 10 derniers
+ for (auto bc : obj.neumanns)
+ out << "\t\t" << *bc << '\n';
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femProblem.h b/cxxfem/src/femProblem.h
new file mode 100644
index 0000000000000000000000000000000000000000..fec3242b5ff48f485a8d3340cb5b25c7343edf8f
--- /dev/null
+++ b/cxxfem/src/femProblem.h
@@ -0,0 +1,50 @@
+#ifndef FEMPROBLEM_H
+#define FEMPROBLEM_H
+
+#include "fem.h"
+#include <map>
+
+namespace fem
+{
+
+/// contains:
+/// - entities: a list of entities
+/// - groups: a list of physical groups
+
+class FEM_API Problem
+{
+public:
+#ifndef SWIG
+ std::vector<Entity *> entities;
+ std::map<std::pair<int, int>, Entity *> entities_by_dimtag;
+
+ std::vector<Group *> groups;
+ std::map<std::pair<int, int>, Group *> groups_by_dimtag;
+ std::map<std::string, Group *> groups_by_name;
+
+ std::vector<std::size_t> nodeTags;
+ std::map<size_t, size_t> nodemap; ///< node tag => index in "nodeTags"
+
+ std::vector<Medium *> media;
+ std::vector<Dirichlet *> dirichlets;
+ std::vector<Neumann *> neumanns;
+
+ std::string quadrature;
+#endif
+public:
+ Problem();
+ ~Problem();
+
+ void sync();
+ void load(std::string const &filename,
+ std::map<std::string, double> const &pars,
+ int order=1);
+
+#ifndef SWIG
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Problem const &obj);
+#endif
+};
+
+} // namespace fem
+
+#endif // FEMPROBLEM_H
diff --git a/cxxfem/src/femSolver.cpp b/cxxfem/src/femSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ec1058bb7f8ffec5816d563a9b0bdaf5737c9d60
--- /dev/null
+++ b/cxxfem/src/femSolver.cpp
@@ -0,0 +1,1008 @@
+#include "femSolver.h"
+#include "femProblem.h"
+#include "femPartition.h"
+#include "femPart.h"
+#include "femDirichlet.h"
+#include "femNeumann.h"
+#include "femGroup.h"
+#include "femEntity.h"
+#include "femElementSet.h"
+#include "femMedium.h"
+#include "femMaps.h"
+#include "OpenMP.h"
+#include "femGmsh.h"
+#include <iomanip>
+#include <algorithm>
+#include <Eigen/IterativeLinearSolvers>
+#ifdef FEM_USE_MKL
+#include <Eigen/PardisoSupport>
+#endif
+
+namespace fem
+{
+
+Solver::Solver(Problem &_pbl) : pbl(_pbl),
+ partition(nullptr),
+ part_Mfree(nullptr),
+ part_Mfixed(nullptr),
+ linear_solver("ldlt")
+{
+ dofs.add("X");
+ dofs.add("Y");
+ dofs.add("Z");
+ bcs.add("FREE");
+ bcs.add("PRESCRIBED");
+
+ std::cout << "#threads = " << OpenMP::get_max_threads() << '\n';
+}
+
+Solver::~Solver()
+{
+ if(debug_level()>0)
+ std::cout << "~Solver()\n";
+ if (partition)
+ delete partition;
+ if (part_Mfree)
+ delete part_Mfree;
+ if (part_Mfixed)
+ delete part_Mfixed;
+
+ // free sparse K
+ for (size_t i = 0; i < K.size(); ++i)
+ for (size_t j = 0; j < K[i].size(); ++j)
+ delete K[i][j];
+
+ // free f, X
+ for (size_t i = 0; i < f.size(); ++i)
+ delete f[i];
+ for (size_t i = 0; i < f.size(); ++i)
+ delete X[i];
+}
+
+/// solve the problem
+
+void Solver::solve()
+{
+ std::cout << pbl;
+
+ create_partition();
+ fill_partition();
+ fill_rhs();
+ fill_lhs();
+ apply_bc();
+ solve_static();
+
+ std::cout << "Summary of CPU times:\n";
+
+ std::ios cout_state(nullptr);
+ cout_state.copyfmt(std::cout);
+
+ std::cout << std::fixed << std::setprecision(1);
+ for (auto const &[name, timer] : timers)
+ {
+ std::cout << '\t' << std::setw(15) << std::left << name << ":"
+ << std::setw(10) << std::right << timer.read() << '\n';
+ }
+
+ std::cout.copyfmt(cout_state);
+}
+
+// create a partition of the dofs
+
+void Solver::create_partition()
+{
+ partition = new Partition();
+ part_Mfree = new Part(*partition, dofs["X"] | dofs["Y"] | dofs["Z"], bcs["FREE"]);
+ part_Mfixed = new Part(*partition, dofs["X"] | dofs["Y"] | dofs["Z"], bcs["PRESCRIBED"]);
+ std::cout << *partition;
+}
+
+// sort nodes & dofs into the defined partition
+// output:
+// dofFixed: is a dof fixed or free?
+// dofPart: in which part is a given dof of a given node?
+// dofIdx: what is its number in its part?
+// partSizes: what are the size of the parts in the partition?
+
+void Solver::fill_partition()
+{
+ timers["fill_partition"].start();
+
+ std::cout << "fill_partition:\n";
+
+ // pbl.nodeTags (size: nb of nodes): gives the gmsh tag for each node
+ // pbl.nodemap (size: max gmsh tag): gives the global node idx (position in nodeTags)
+ // from the gmsh tag
+ // (some space could be unused if non-contiguous numbering)
+ //
+ // dofPart (size: nb of nodes x size of dofs): for each node/dof, gives its partition number
+ // dofIdx (size: nb of nodes x size of dofs): for each node/dof, gives its index in its partition
+
+ // fill the 'dofFixed' array:
+ // T X Y Z
+ // node 0 [ x x x x ]
+ // node 1 [ x x x x ] where 'x' is PRESCRIBED or FREE
+ // node 2 [ x x x x ]
+ // ...
+
+ dofFixed.setConstant(pbl.nodeTags.size(), dofs.size(), 0); // = false == FREE by default
+
+ for (auto &bc : pbl.dirichlets)
+ {
+ DOF bc_dof = dofs[bc->m_dofname];
+ for (auto &entity : bc->group->entities)
+ for (auto &elems : entity->esets)
+ for (size_t i = 0; i < elems->etags.size(); ++i)
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ // size_t n = elems->ntags_v[i * elems->type_nnods + j];
+ size_t n = elems->ntags(i, j);
+ size_t gidx = pbl.nodemap[n];
+ dofFixed(gidx, bc_dof.index()) = 1; // PRESCRIBED
+ }
+ }
+
+ Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "dofFixed =\n"
+ // << dofFixed.format(fmt) << '\n';
+
+ // fill dofPart (default:-1 => no partition assigned)
+ // T X Y Z
+ // node 0 [ x x x x ]
+ // node 1 [ x x x x ] where 'x' is the partition of the dof
+ // node 2 [ x x x x ]
+ // ...
+
+ dofPart.setConstant(pbl.nodeTags.size(), dofs.size(), -1);
+
+ for (size_t i = 0; i < pbl.nodeTags.size(); ++i)
+ for (size_t j = 0; j < dofs.size(); ++j) // TODO: iterator
+ {
+ DOF dof = dofs[j];
+ DOF bc = dofFixed(i, j) ? bcs["PRESCRIBED"] : bcs["FREE"];
+ for (auto &p : partition->parts)
+ if ((p->dofmask & dof) && (p->bcmask & bc))
+ {
+ dofPart(i, j) = p->idx;
+ break;
+ }
+ }
+
+ // std::cout << "dofPart =\n"
+ // << dofPart.format(fmt) << '\n';
+
+ // count nodes in partitions
+ int nnods_Mfree = (dofPart == part_Mfree->idx).count();
+ std::cout << "\tfree M dofs = " << nnods_Mfree << '\n';
+ int nnods_Mfixed = (dofPart == part_Mfixed->idx).count();
+ std::cout << "\tfixed M dofs = " << nnods_Mfixed << std::endl;
+
+ // renumber nodes in each partition
+ // fill dofIdx (default:-1 => no part assigned)
+ // T X Y Z
+ // node 0 [ x x x x ]
+ // node 1 [ x x x x ] where 'x' is the index of the dof in its partition
+ // node 2 [ x x x x ]
+ // ...
+ // fill also 'partSizes' (nb of dofs in each part)
+
+ dofIdx.setConstant(pbl.nodeTags.size(), dofs.size(), 0);
+ partSizes.resize(partition->parts.size(), 0);
+
+ for (size_t i = 0; i < pbl.nodeTags.size(); ++i)
+ for (size_t j = 0; j < dofs.size(); ++j) // TODO: iterator
+ {
+ int p = dofPart(i, j);
+ dofIdx(i, j) = partSizes[p];
+ partSizes[p]++;
+ }
+
+ // std::cout << "dofIdx =\n"
+ // << dofIdx.format(fmt) << '\n';
+ // std::cout << "partSizes = " << partSizes << std::endl;
+
+ // TODO: check that all dofs are assigned?
+ timers["fill_partition"].stop();
+ std::cout << "\tcpu = " << timers["fill_partition"].read() << "s\n";
+}
+
+/// fill right-hand side of the system of equations
+
+void Solver::fill_rhs()
+{
+ timers["fill_rhs"].start();
+
+ std::cout << "assembly of LHS...\n";
+
+ // -- prepare structural matrices
+ //
+ // Kijv[nparts][nparts], Cijv[nparts][nparts]:
+ // "ijv" versions of K, C and M used during assembly
+
+ std::vector<std::vector<std::vector<Eigen::Triplet<double>> *>> Kijv(partition->parts.size());
+ for (size_t i = 0; i < partition->parts.size(); ++i)
+ {
+ Kijv[i].resize(partition->parts.size());
+ for (size_t j = 0; j < partition->parts.size(); ++j)
+ Kijv[i][j] = new std::vector<Eigen::Triplet<double>>;
+ }
+
+ // -- pre-compute sub-matrices structures
+ std::vector<std::vector<size_t>> ntriplets(partition->parts.size(), std::vector<size_t>(partition->parts.size(), 0));
+
+ std::vector<DOF> dofs_M = {dofs["X"], dofs["Y"], dofs["Z"]};
+ for (auto medium : pbl.media)
+ for (auto entity : medium->group->entities)
+ for (auto elems : entity->esets)
+ {
+ const int ndim = elems->type_dim;
+ for (size_t ie = 0; ie < elems->etags.size(); ++ie)
+ {
+ // brutal assembly (needs optimisation)
+ for (int i = 0; i < elems->type_nnods; ++i)
+ {
+ int i_tag = elems->ntags(ie, i); // gmsh tag of node i
+ int i_gidx = pbl.nodemap[i_tag]; // global index of node i
+ auto i_part = dofPart.row(i_gidx); // parts of dofs of node i
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ int j_tag = elems->ntags(ie, j); // gmsh tag of node j
+ int j_gidx = pbl.nodemap[j_tag]; // global index of node j
+ auto j_part = dofPart.row(j_gidx); // parts of dofs of node j
+
+ for (int ni = 0; ni < ndim; ++ni)
+ {
+ DOF &di = dofs_M[ni];
+ for (int nj = 0; nj < ndim; ++nj)
+ {
+ DOF &dj = dofs_M[nj];
+
+ int ip = i_part(di.index());
+ int jp = j_part(dj.index());
+ ntriplets[ip][jp]++;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // reserve the correct amount of triplets
+ // (avoids reallocation in later push_backs)
+ for (size_t i = 0; i < partition->parts.size(); ++i)
+ {
+ Kijv[i].resize(partition->parts.size());
+ for (size_t j = 0; j < partition->parts.size(); ++j)
+ {
+ std::cout << "expecting " << ntriplets[i][j]
+ << " triplets for submatrix " << i << "," << j << '\n';
+ Kijv[i][j]->reserve(ntriplets[i][j]);
+ }
+ }
+
+ // -- compute elemental matrices
+
+ for (auto medium : pbl.media)
+ for (auto entity : medium->group->entities)
+ for (auto elems : entity->esets)
+ {
+ std::cout << "\tprocessing " << *elems << '\n';
+ size_t nelems = elems->etags.size();
+ std::cout << "\t\tnelems = " << nelems << '\n';
+ // integration points
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(elems->type_id,
+ pbl.quadrature,
+ localCoord,
+ weights);
+ size_t nGP = weights.size();
+ std::cout << "\t\tnGP = " << nGP << '\n';
+ // jacobians
+ std::vector<double> jacobians_v;
+ std::vector<double> determinants_v;
+ std::vector<double> coord_v;
+ gmsh::model::mesh::getJacobians(elems->type_id, localCoord,
+ jacobians_v, determinants_v,
+ coord_v, entity->tag);
+ // reshape (map) into matrices
+ JacoMap jacobians(jacobians_v, nelems, nGP);
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+ determinants(&determinants_v[0], nelems, nGP);
+ CoordMap coord(coord_v, nelems, nGP);
+
+ // print GP data
+ // for (size_t i = 0; i < nelems; ++i)
+ // {
+ // std::cout << "\telem #" << i << " (tag=" << elems->ntags[i] << "):\n";
+ // for (size_t j = 0; j < nGP; ++j)
+ // {
+ // std::cout << "\t\tGP #" << j + 1 << '/' << nGP << ":\n";
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t\t\t\t[", "]");
+ // std::cout << "\t\t\tJ =\n"
+ // << jacobians(i, j).format(fmt) << '\n';
+ // std::cout << "\t\t\tdet = " << jacobians(i, j).determinant()
+ // << " (gmsh = " << determinants(i, j) << ")\n";
+ // std::cout << "\t\t\tcoord =\n"
+ // << coord(i, j).format(fmt) << '\n';
+ // }
+ // }
+
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "determinants =\n"
+ // << determinants.format(fmt) << '\n';
+ // std::cout << "determinants_v =\n" << determinants_v << '\n';
+ // std::cout << "coord_v =\n" << coord_v << '\n';
+ // std::cout << "coord_v.size() =" << coord_v.size() << '\n';
+
+ // basis functions
+ int numComponents; // unused
+ std::vector<double> basisF_v;
+ int numOrientations; // unused
+ gmsh::model::mesh::getBasisFunctions(elems->type_id,
+ localCoord, "Lagrange",
+ numComponents, // unused
+ basisF_v,
+ numOrientations); // unused
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+ basisF(&basisF_v[0], nGP, elems->type_nnods);
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "basisF =\n"
+ // << basisF.format(fmt) << '\n';
+ // std::cout << "basisF_v =\n"
+ // << basisF_v << '\n';
+ // exit(1);
+
+ // basis functions (derivatives)
+ std::vector<double> basisDF_v;
+ gmsh::model::mesh::getBasisFunctions(elems->type_id,
+ localCoord, "GradLagrange",
+ numComponents, // reused
+ basisDF_v,
+ numOrientations); // reused
+ DShapeMap basisDF(basisDF_v, nGP, elems->type_nnods);
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "basisDF =\n"
+ // << basisDF(0).format(fmt) << '\n';
+ // std::cout << "basisDF_v =\n"
+ // << basisDF_v << '\n';
+ // exit(1);
+
+ #ifdef _MSC_VER
+ #define OPENMP_INT int
+ #else
+ #define OPENMP_INT size_t
+ #endif
+
+ // loop over the elements
+ #pragma omp parallel for
+ for (OPENMP_INT ie = 0; ie < elems->etags.size(); ++ie)
+ {
+ // std::cout << "\tprocessing element #" << elems->etags[ie] << '\n';
+ const int ndim = elems->type_dim;
+
+ // compute elemental matrices
+ Eigen::MatrixXd Ke_MM = Eigen::MatrixXd::Zero(ndim * elems->type_nnods, ndim * elems->type_nnods);
+
+ // loop over Gauss points
+ for (size_t igp = 0; igp < nGP; ++igp)
+ {
+ auto N_T = basisF.row(igp); // matrix of shape functions; shape=(1,elems.type_nnods)
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "N_T (igp=" << igp << ") = " << N_T.format(fmt) << '\n';
+ auto DN = basisDF(igp).block(0, 0, elems->type_nnods, ndim); // matrix of shp fct derivatives in 'xi' space; shape=(elems.type_nnods,3)
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "DN (igp=" << igp << ") =\n"
+ // << DN.format(fmt) << '\n';
+ auto J = jacobians(ie, igp).block(0, 0, ndim, ndim);
+ auto Jt = J.transpose(); // transposed jacobian matrix; shape=(3,3)
+ auto Jtinv = Jt.inverse(); // transposed jacobian inverse; shape=(3,3)
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "Jt*Jtinv = " << (Jt*Jtinv).format(fmt) << '\n';
+ auto B_T = Jtinv * DN.transpose(); // matrix of shp fct derivatives in 'x' space; shape=(3, nnods_by_element)
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "B_T = " << B_T.format(fmt) << '\n';
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "basisF = " << basisF.format(fmt) << '\n';
+
+ Eigen::MatrixXd N_M(ndim, elems->type_nnods * ndim);
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ // double v = basisF(igp, j);
+ double v = N_T(j);
+ // std::cout << "v = " << v << '\n';
+ N_M.block(0, ndim * j, ndim, ndim) = Eigen::MatrixXd::Identity(ndim, ndim) * v;
+ }
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "N_M = " << N_M.format(fmt) << '\n';
+
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "B_T = " << B_T.format(fmt) << '\n';
+
+ Eigen::MatrixXd B_M;
+
+ switch (ndim)
+ {
+ case 2:
+
+ B_M.resize(3, elems->type_nnods * ndim); // 3 en 2D; 6 en 3D
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ double ddx = B_T(0, j);
+ double ddy = B_T(1, j);
+ Eigen::MatrixXd Bs(3, 2);
+ Bs << ddx, 0.0,
+ 0.0, ddy,
+ ddy, ddx;
+ B_M.block(0, ndim * j, 3, 2) = Bs;
+ }
+ break;
+ case 3:
+
+ B_M.resize(6, elems->type_nnods * ndim); // 3 en 2D; 6 en 3D
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ // double v = basisF(igp, j);
+ double ddx = B_T(0, j);
+ double ddy = B_T(1, j);
+ double ddz = B_T(2, j);
+ Eigen::MatrixXd Bs(6, 3);
+ Bs << ddx, 0.0, 0.0,
+ 0.0, ddy, 0.0,
+ 0.0, 0.0, ddz,
+ ddy, ddx, 0.0,
+ ddz, 0.0, ddx,
+ 0.0, ddz, ddy;
+ // std::cout << "v = " << v << '\n';
+ B_M.block(0, ndim * j, 6, 3) = Bs;
+ }
+ break;
+ default:
+ throw std::runtime_error("dimension " + std::to_string(ndim) + " not implemented");
+ }
+ // std::cout << "B_M = " << B_M.format(fmt) << '\n';
+
+ Ke_MM += B_M.transpose() * medium->hooke(ndim) * B_M * determinants(ie, igp) * weights[igp];
+ //
+ // exit(1);
+ }
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "Ke_MM = \n"
+ // << Ke_MM.format(fmt) << '\n';
+ // exit(1);
+
+ std::vector<DOF> dofs_M = {dofs["X"], dofs["Y"], dofs["Z"]};
+
+ // brutal assembly (needs optimisation)
+ #pragma omp critical
+ for (int i = 0; i < elems->type_nnods; ++i)
+ {
+ int i_tag = elems->ntags(ie, i); // gmsh tag of node i
+ int i_gidx = pbl.nodemap[i_tag]; // global index of node i
+ auto i_part = dofPart.row(i_gidx); // parts of dofs of node i
+ auto i_idx = dofIdx.row(i_gidx); // index of these dofs in their parts
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ int j_tag = elems->ntags(ie, j); // gmsh tag of node j
+ int j_gidx = pbl.nodemap[j_tag]; // global index of node j
+ auto j_part = dofPart.row(j_gidx); // parts of dofs of node j
+ auto j_idx = dofIdx.row(j_gidx); // index of these dofs in their parts
+
+ for (int ni = 0; ni < ndim; ++ni)
+ {
+ DOF &di = dofs_M[ni];
+ for (int nj = 0; nj < ndim; ++nj)
+ {
+ DOF &dj = dofs_M[nj];
+ // std::cout << "(" << di.index() << ',' << dj.index() << ")\n";
+ // get vector of triplets
+ auto K_sp = Kijv[i_part(di.index())][j_part(dj.index())];
+ K_sp->push_back(Eigen::Triplet<double>(
+ i_idx(di.index()),
+ j_idx(dj.index()),
+ Ke_MM(i * ndim + ni, j * ndim + nj)));
+ }
+ }
+ }
+ } // assembly
+
+ } // loop over the elements
+ } // for (auto &elems : entity->esets)
+
+ // create sparse K
+
+ K.resize(partition->parts.size());
+ for (size_t i = 0; i < partition->parts.size(); ++i)
+ {
+ K[i].resize(partition->parts.size());
+ for (size_t j = 0; j < partition->parts.size(); ++j)
+ {
+ K[i][j] = new Eigen::SparseMatrix<double>(partSizes[i], partSizes[j]);
+ K[i][j]->setFromTriplets(Kijv[i][j]->begin(), Kijv[i][j]->end()); // TODO: use std::list!
+ }
+ }
+
+ // free triplets
+ for (size_t i = 0; i < Kijv.size(); ++i)
+ for (size_t j = 0; j < Kijv[i].size(); ++j)
+ delete Kijv[i][j];
+
+ timers["fill_rhs"].stop();
+ std::cout << "\tcpu = " << timers["fill_rhs"].read() << "s\n";
+}
+
+/// compute right-hand side 'f' (sources)
+
+void Solver::fill_lhs()
+{
+ timers["fill_lhs"].start();
+
+ std::cout << "assembly of RHS...\n";
+
+ // create RHS
+ f.resize(partition->parts.size());
+ for (size_t i = 0; i < partition->parts.size(); ++i)
+ {
+ f[i] = new Eigen::VectorXd(partSizes[i]);
+ f[i]->setZero();
+ }
+
+ // sources
+ // ...
+
+ // std::cout << "Neumann BCs...\n";
+
+ for (auto src : pbl.neumanns)
+ {
+ DOF bc_dof = dofs[src->m_dofname];
+ if (src->group)
+ {
+ for (auto entity : src->group->entities) // each source has a target group of entities
+ for (auto elems : entity->esets) // each entity has a several sets of elements
+ {
+ size_t nelems = elems->etags.size();
+ // integration points
+ std::vector<double> localCoord;
+ std::vector<double> weights;
+ gmsh::model::mesh::getIntegrationPoints(elems->type_id,
+ pbl.quadrature,
+ localCoord,
+ weights);
+ size_t nGP = weights.size();
+ // jacobians
+ std::vector<double> jacobians_v;
+ std::vector<double> determinants_v;
+ std::vector<double> coord_v;
+ gmsh::model::mesh::getJacobians(elems->type_id, localCoord,
+ jacobians_v, determinants_v,
+ coord_v, entity->tag);
+ // reshape (map) into matrices
+ JacoMap jacobians(jacobians_v, nelems, nGP);
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+ determinants(&determinants_v[0], nelems, nGP);
+ CoordMap coord(coord_v, nelems, nGP);
+ // basis functions
+ int numComponents; // unused
+ std::vector<double> basisF_v;
+ int numOrientations; // unused
+ gmsh::model::mesh::getBasisFunctions(elems->type_id,
+ localCoord, "Lagrange",
+ numComponents, // unused
+ basisF_v,
+ numOrientations); // unused
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+ basisF(&basisF_v[0], nGP, elems->type_nnods);
+
+ // loop over the elements
+ for (size_t ie = 0; ie < elems->etags.size(); ++ie)
+ {
+ Eigen::VectorXd fe = Eigen::VectorXd::Zero(elems->type_nnods);
+ // loop over Gauss points
+ for (size_t igp = 0; igp < nGP; ++igp)
+ {
+ auto N_T = basisF.row(igp); // matrix of shape functions; shape=(1,elems.type_nnods)
+ fe += src->m_val * N_T * determinants(ie, igp) * weights[igp];
+ }
+
+ // brutal assembly (needs optimisation)
+ for (int i = 0; i < elems->type_nnods; ++i)
+ {
+ int i_tag = elems->ntags(ie, i); // gmsh tag of node i
+ int i_gidx = pbl.nodemap[i_tag]; // global index of node i
+ auto i_part = dofPart.row(i_gidx)[bc_dof.index()]; // parts of dofs of node i
+ auto i_idx = dofIdx.row(i_gidx)[bc_dof.index()]; // index of these dofs in their parts
+ (*f[i_part])[i_idx] += fe[i];
+ }
+ }
+ }
+ } // if (src->group)
+ else
+ {
+ // neumann applied on a single node
+ int i_tag = src->m_nodetag; // gmsh tag of node i
+ int i_gidx = pbl.nodemap[i_tag]; // global index of node i
+ auto i_part = dofPart.row(i_gidx)[bc_dof.index()]; // parts of dofs of node i
+ auto i_idx = dofIdx.row(i_gidx)[bc_dof.index()]; // index of these dofs in their parts
+ (*f[i_part])[i_idx] += src->m_val;
+ }
+ }
+
+ timers["fill_lhs"].stop();
+ std::cout << "\tcpu = " << timers["fill_lhs"].read() << "s\n";
+
+}
+
+// create the X vector (containing dofs)
+// and prescribe dirichlet boundary conditions
+void Solver::apply_bc()
+{
+ timers["apply_bc"].start();
+ std::cout << "applying dirichlet boundary conditions...\n";
+
+ // create vectors of dofs
+ X.resize(partition->parts.size());
+ for (size_t i = 0; i < partition->parts.size(); ++i)
+ {
+ X[i] = new Eigen::VectorXd(partSizes[i]);
+ X[i]->setZero();
+ }
+
+ // apply boundary conditions
+
+ for (auto &bc : pbl.dirichlets)
+ {
+ DOF bc_dof = dofs[bc->m_dofname];
+ for (auto &entity : bc->group->entities)
+ for (auto &elems : entity->esets)
+ for (size_t i = 0; i < elems->etags.size(); ++i)
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ // size_t n = elems->ntags_v[i * elems->type_nnods + j];
+ size_t n = elems->ntags(i, j);
+ size_t gidx = pbl.nodemap[n];
+ int part = dofPart(gidx, bc_dof.index());
+ int idx = dofIdx(gidx, bc_dof.index());
+ (*X[part])[idx] = bc->m_val;
+ }
+ }
+
+ timers["apply_bc"].stop();
+ std::cout << "\tcpu = " << timers["apply_bc"].read() << "s\n";
+}
+
+void Solver::solve_static()
+{
+ timers["solve_static"].start();
+
+ std::cout << "solving linear system...\n";
+
+ // create vectors of tmp_rhs
+ std::vector<Eigen::VectorXd *> tmp_rhs(partition->parts.size());
+ for (size_t i = 0; i < partition->parts.size(); ++i)
+ {
+ tmp_rhs[i] = new Eigen::VectorXd(partSizes[i]);
+ tmp_rhs[i]->setZero();
+ }
+
+ // system: (F="free" ; P="prescribed")
+ // [ K_FF K_FP ] [ x_F ] = [ f_F ]
+ // [ K_PF K_PP ] [ x_P ] [ f_P ]
+ // construct right-hand sides
+ // [ K_FF ][ x_F ] + [ K_FP ][ x_P ] = [ f_F ]
+ // => [ K_FF ][ x_F ] = [ f_F ] - [ K_FP ][ x_P ]
+ // <-------tmp_rhs_F------->
+ // (same procedure for both thermal and mechanical problems)
+
+ size_t MF = part_Mfree->idx;
+ size_t MP = part_Mfixed->idx;
+
+ // mechanical
+ *tmp_rhs[MF] = *f[MF] - (*K[MF][MP]) * (*X[MP]);
+
+ // mechanical problem: [ x_F ] = [ K_FF ]^(-1) [ tmp_rhs_F ]
+ std::cout << "\t#unknowns = " << K[MF][MF]->rows() << '\n';
+ std::cout << "\t#threads = " << OpenMP::get_max_threads() << '\n';
+
+ if (linear_solver == "cg")
+ {
+ // --conjugate gradient
+ std::cout << "\tusing Eigen::ConjugateGradient (linear_solver=\""+linear_solver+"\")\n";
+ Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower | Eigen::Upper> solver;
+ solver.compute(*K[MF][MF]);
+ if(solver.info()!=Eigen::Success)
+ throw std::runtime_error("solver.compute() returned "+std::to_string(solver.info()));
+ *X[MF] = solver.solve(*tmp_rhs[MF]);
+ std::cout << "\t#iterations: " << solver.iterations() << std::endl;
+ std::cout << "\testimated error: " << solver.error() << std::endl;
+ if(solver.info()!=Eigen::Success)
+ throw std::runtime_error("solver.solve() returned "+std::to_string(solver.info()));
+ }
+ else if (linear_solver == "ldlt")
+ {
+ // --direct solver (LDLT)
+ std::cout << "\tusing Eigen::SimplicialLDLT (linear_solver=\""+linear_solver+"\")\n";
+ Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+ solver.compute(*K[MF][MF]);
+ if(solver.info()!=Eigen::Success)
+ throw std::runtime_error("solver.compute() returned "+std::to_string(solver.info()));
+ *X[MF] = solver.solve(*tmp_rhs[MF]); // use the factorization to solve for the given right-hand side
+ if(solver.info()!=Eigen::Success)
+ throw std::runtime_error("solver.solve() returned "+std::to_string(solver.info()));
+ }
+#ifdef FEM_USE_MKL
+ else if (linear_solver == "pardiso")
+ {
+ // -- symmetric pardiso (LDLT)
+ std::cout << "\tusing Eigen::PardisoLDLT (linear_solver=\""+linear_solver+"\")\n";
+ Eigen::PardisoLDLT<Eigen::SparseMatrix<double>> solver;
+ solver.compute(*K[MF][MF]);
+ if(solver.info()!=Eigen::Success)
+ throw std::runtime_error("solver.compute() returned "+std::to_string(solver.info()));
+ *X[MF] = solver.solve(*tmp_rhs[MF]); // use the factorization to solve for the given right-hand side
+ if(solver.info()!=Eigen::Success)
+ throw std::runtime_error("solver.solve() returned "+std::to_string(solver.info()));
+ }
+#endif // FEM_USE_MKL
+ else
+ {
+ throw std::runtime_error("unknown linear_solver: "+linear_solver);
+ }
+
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "X[MF] = \n"
+ // << (*X[MF]).format(fmt) << '\n';
+
+ // compute reactions: [ f_P ] = [ K_PF ][ x_F ] + [ K_PP ][ x_P ]
+ // values written previously in [ f_P ] during source / Neuman-BCs
+ // treatment are overwritten
+
+ *f[MP] = (*K[MP][MF]) * (*X[MF]) + (*K[MP][MP]) * (*X[MP]);
+
+ // std::cout << "f[MP] = \n"
+ // << (*f[MP]).format(fmt) << '\n';
+
+ // energy: work of ext forces (using their "final" value - thus twice the int energy)
+ double workF = (*f[MP]).dot(*X[MP]) + (*f[MF]).dot(*X[MF]);
+ std::cout << "workF = " << workF << '\n';
+ int_energy = workF/2;
+
+ // free vectors of tmp_rhs
+ for (size_t i = 0; i < tmp_rhs.size(); ++i)
+ delete tmp_rhs[i];
+
+ timers["solve_static"].stop();
+ std::cout << "\tcpu = " << timers["solve_static"].read() << "s\n";
+}
+
+// strain and stresses (epsilon, sigma) per node, per element
+// (discontinuous across element boundaries)
+
+void Solver::compute_sig_epl(std::vector<size_t> &etags_all, // tags of elements
+ std::vector<std::vector<double>> &sigma_all, // values of stresses at the nodes of each element
+ std::vector<std::vector<double>> &epsilon_all, // values of strains at the nodes of each element
+ std::map<size_t, std::vector<Eigen::VectorXd>> &sig_by_node, // value of all the stress vectors (Voigt) at a given node (key of the dict)
+ std::map<size_t, std::vector<Eigen::VectorXd>> &eps_by_node) // same for strain vector
+{
+ timers["compute_sig_epl"].start();
+
+ std::cout << "computing strain and stresses...\n";
+
+ // pre-allocate things
+ std::map<size_t, size_t> nneighbours;
+ size_t nelems=0;
+ for (auto &medium : pbl.media)
+ for (auto &entity : medium->group->entities)
+ {
+ for (auto &elems : entity->esets)
+ {
+ nelems+=elems->etags.size(); // accumulate nb of elems
+
+ for (size_t ie = 0; ie < elems->etags.size(); ++ie)
+ for (int inod = 0; inod < elems->ntags.cols(); ++inod)
+ {
+ size_t ntag = elems->ntags(ie, inod);
+ nneighbours[ntag]++;
+ }
+ }
+ }
+ std::cout << "\texpecting " << nelems << " elements\n";
+ etags_all.reserve(nelems);
+ epsilon_all.reserve(nelems);
+ sigma_all.reserve(nelems);
+
+ size_t maxneigh=0;
+ for(auto const &[ntag, nneigh] : nneighbours)
+ {
+ maxneigh = std::max(nneigh, maxneigh);
+ eps_by_node[ntag].reserve(nneigh);
+ sig_by_node[ntag].reserve(nneigh);
+ }
+ std::cout << "\tmax #neighbours = " << maxneigh << '\n';
+
+ // loop over the media
+
+ for (auto &medium : pbl.media)
+ for (auto &entity : medium->group->entities)
+ for (auto &elems : entity->esets)
+ {
+ // std::cout << "processing " << *elems << '\n';
+ size_t nelems = elems->etags.size();
+
+ // jacobians at nodes
+ std::vector<double> jacobians_v;
+ std::vector<double> determinants_v;
+ std::vector<double> coord_v;
+ gmsh::model::mesh::getJacobians(elems->type_id, elems->type_localNodeCoord,
+ jacobians_v, determinants_v,
+ coord_v, entity->tag);
+ // reshape (map) into matrices
+ JacoMap jacobians(jacobians_v, nelems, elems->type_nnods);
+ Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
+ determinants(&determinants_v[0], nelems, elems->type_nnods);
+ CoordMap coord(coord_v, nelems, elems->type_nnods);
+
+ // basis functions (derivatives)
+ int numComponents; // unused
+ int numOrientations; // unused
+ std::vector<double> basisDF_v;
+ gmsh::model::mesh::getBasisFunctions(elems->type_id,
+ elems->type_localNodeCoord,
+ "GradLagrange",
+ numComponents, // unused
+ basisDF_v,
+ numOrientations); // unused
+ DShapeMap basisDF(basisDF_v, elems->type_nnods, elems->type_nnods);
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "basisDF =\n"
+ // << basisDF(0).format(fmt) << '\n';
+ // std::cout << "basisDF_v =\n"
+ // << basisDF_v << '\n';
+ // exit(1);
+
+ // loop over the elements
+ //#pragma omp parallel
+ for (size_t ie = 0; ie < elems->etags.size(); ++ie)
+ {
+ // std::cout << "element " << ie+1 << "/" << elems->etags.size() << '\n';
+ size_t etag = elems->etags[ie];
+ const int ndim = elems->type_dim;
+
+ // compute the nodal displacement vectors ("displ")
+ // and the temperatures ("temp")
+ // of the current element
+ Eigen::VectorXd displ(elems->ntags.cols() * ndim);
+ // loop over the nodes
+ for (int inod = 0; inod < elems->ntags.cols(); ++inod)
+ {
+ size_t ntag = elems->ntags(ie, inod);
+ size_t gidx = pbl.nodemap[ntag];
+ for (int di = 0; di < ndim; ++di)
+ {
+ DOF dof = dofs[di];
+ int part = dofPart(gidx, dof.index());
+ int idx = dofIdx(gidx, dof.index());
+ displ[ndim * inod + di] = (*X[part])[idx];
+ }
+ }
+
+ std::vector<double> eps_nodes(9 * elems->ntags.cols(), 0.0); // strains at the nodes of the element in a 1d array
+ std::vector<double> sig_nodes(9 * elems->ntags.cols(), 0.0); // stresses at the nodes of the element in a 1d array
+
+ // loop over the nodes
+ // and compute the stresses/strains/heat_fluxes at the nodes
+
+ for (int inod = 0; inod < elems->ntags.cols(); ++inod)
+ {
+ size_t ntag = elems->ntags(ie, inod);
+ auto DN = basisDF(inod).block(0, 0, elems->type_nnods, ndim); // matrix of shp fct derivatives in 'xi' space; shape=(elems.type_nnods,3)
+ auto J = jacobians(ie, inod).block(0, 0, ndim, ndim);
+ auto Jt = J.transpose(); // transposed jacobian matrix; shape=(3,3)
+ auto Jtinv = Jt.inverse(); // transposed jacobian inverse; shape=(3,3)
+ auto B_T = Jtinv * DN.transpose(); // matrix of shp fct derivatives in 'x' space; shape=(3, nnods_by_element)
+
+ Eigen::MatrixXd B_M;
+ switch (ndim)
+ {
+ case 2:
+ B_M.resize(3, elems->type_nnods * ndim); // 3 en 2D; 6 en 3D
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ double ddx = B_T(0, j);
+ double ddy = B_T(1, j);
+ Eigen::MatrixXd Bs(3, 2);
+ Bs << ddx, 0.0,
+ 0.0, ddy,
+ ddy, ddx;
+ B_M.block(0, ndim * j, 3, 2) = Bs;
+ }
+ break;
+ case 3:
+ B_M.resize(6, elems->type_nnods * ndim); // 3 en 2D; 6 en 3D
+ for (int j = 0; j < elems->type_nnods; ++j)
+ {
+ // double v = basisF(igp, j);
+ double ddx = B_T(0, j);
+ double ddy = B_T(1, j);
+ double ddz = B_T(2, j);
+ Eigen::MatrixXd Bs(6, 3);
+ Bs << ddx, 0.0, 0.0,
+ 0.0, ddy, 0.0,
+ 0.0, 0.0, ddz,
+ ddy, ddx, 0.0,
+ ddz, 0.0, ddx,
+ 0.0, ddz, ddy;
+ // std::cout << "v = " << v << '\n';
+ B_M.block(0, ndim * j, 6, 3) = Bs;
+ }
+ break;
+ default:
+ throw std::runtime_error("dimension " + std::to_string(ndim) + " not implemented");
+ }
+
+ auto eps = B_M * displ; // nodal epsilon (Voigt)
+ auto sig = medium->hooke(ndim) * eps; // nodal sigma (Voigt)
+ // Eigen::IOFormat fmt(6, 0, ", ", "\n", "\t[", "]");
+ // std::cout << "sig = \n"
+ // << sig.format(fmt) << '\n';
+
+ // convert to 3D Voigt vectors
+ Eigen::VectorXd eps3d(6), sig3d(6); // TODO: use fixed-size vectors
+ if (ndim == 3)
+ {
+ eps3d = eps;
+ sig3d = sig;
+ }
+ else
+ {
+ double eps_zz = -medium->m_nu / medium->m_E * (sig[0] + sig[1]);
+ eps3d << eps(0), eps(1), eps_zz,
+ eps(2), 0.0, 0.0;
+ sig3d << sig(0), sig(1), 0.0,
+ sig(2), 0.0, 0.0;
+ }
+
+ // [ [0] [ [XX] [ [0] [3] [4] [ [XX] [XY] [XZ]
+ // [1] [YY] [3] [1] [5] [XY] [YY] [YZ]
+ // [2] [ZZ] [4] [5] [2] ] [XZ] [YZ] [ZZ] ]
+ // [3] [XY]
+ // [4] [XZ]
+ // [5] ] [YZ] ]
+
+ std::vector<int> voigt3d = {0, 3, 4, 3, 1, 5, 4, 5, 2};
+ for (int i = 0; i < 9; ++i)
+ {
+ eps_nodes.at(inod * 9 + i) = eps3d(voigt3d[i]); // TODO: remove .at => []
+ sig_nodes.at(inod * 9 + i) = sig3d(voigt3d[i]);
+ }
+
+ //#pragma omp critical
+ {
+ // store data also by node for smoothed fields
+ eps_by_node[ntag].push_back(eps3d);
+ sig_by_node[ntag].push_back(sig3d);
+ }
+
+ } // for (int inod = 0; ... (nodes of element ie)
+
+ //#pragma omp critical
+ {
+ etags_all.push_back(etag);
+ epsilon_all.push_back(eps_nodes);
+ sigma_all.push_back(sig_nodes);
+ }
+
+ } // for (size_t ie = 0; ... (elements)
+ } // for (auto &elems : entity->esets)
+
+ timers["compute_sig_epl"].stop();
+ std::cout << "\tcpu = " << timers["compute_sig_epl"].read() << "s\n";
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Solver const &obj)
+{
+ out << "Solver\n";
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femSolver.h b/cxxfem/src/femSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..52af8f36f295e5559899745d68e2afe3fa05fc64
--- /dev/null
+++ b/cxxfem/src/femSolver.h
@@ -0,0 +1,69 @@
+#ifndef FEMSOLVER_H
+#define FEMSOLVER_H
+
+#include "fem.h"
+#include "femDOF.h"
+#include "femTimer.h"
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace fem
+{
+
+/// FEM solver
+
+class FEM_API Solver
+{
+public:
+#ifndef SWIG
+ Problem &pbl; ///< problem data
+ Partition *partition; ///< partition of the dofs
+ Part *part_Mfree;
+ Part *part_Mfixed;
+
+ DOFs dofs;
+ DOFs bcs;
+
+ Eigen::ArrayXXi dofFixed; ///< is a dof fixed or free?
+ Eigen::ArrayXXi dofPart; ///< in which part is a given dof of a given node?
+ Eigen::ArrayXXi dofIdx;
+ std::vector<size_t> partSizes;
+
+ std::vector<std::vector<Eigen::SparseMatrix<double> *>> K;
+ std::vector<Eigen::VectorXd *> f;
+ std::vector<Eigen::VectorXd *> X;
+
+ std::map<std::string, Timer> timers;
+
+#endif
+ std::string linear_solver;
+ double int_energy;
+
+public:
+ Solver(Problem &_pbl);
+ ~Solver();
+
+ void solve();
+
+#ifndef SWIG
+ void compute_sig_epl(std::vector<size_t> &etags_all,
+ std::vector<std::vector<double>> &sigma_all,
+ std::vector<std::vector<double>> &epsilon_all,
+ std::map<size_t, std::vector<Eigen::VectorXd>> &sig_by_node,
+ std::map<size_t, std::vector<Eigen::VectorXd>> &eps_by_node);
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Solver const &obj);
+#endif
+
+private:
+ void create_partition();
+ void fill_partition();
+ void fill_rhs();
+ void fill_lhs();
+ void apply_bc();
+ void solve_static();
+};
+
+} // namespace fem
+
+#endif // FEMSOLVER_H
diff --git a/cxxfem/src/femTimer.cpp b/cxxfem/src/femTimer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bcda2e44b4bf5b51835f999777b0277ec00e0d01
--- /dev/null
+++ b/cxxfem/src/femTimer.cpp
@@ -0,0 +1,55 @@
+#include "femTimer.h"
+
+namespace fem
+{
+
+Timer::Timer()
+{
+ reset();
+}
+
+Timer::~Timer()
+{
+ // std::cout << "~Timer()\n";
+}
+
+void Timer::start()
+{
+ running = true;
+ t0 = std::chrono::high_resolution_clock::now();
+}
+
+void Timer::stop()
+{
+ telapsed = read();
+ t0 = std::chrono::high_resolution_clock::time_point();
+ running = false;
+}
+
+void Timer::reset()
+{
+ t0 = std::chrono::high_resolution_clock::time_point();
+ telapsed = 0.0;
+ running = false;
+}
+
+double Timer::read() const
+{
+ if (running)
+ {
+ auto tnow = std::chrono::high_resolution_clock::now();
+ double dt = std::chrono::duration_cast<std::chrono::duration<double>>(tnow - t0).count();
+ return telapsed + dt;
+ }
+ else
+ return telapsed;
+}
+
+FEM_API std::ostream &
+operator<<(std::ostream &out, Timer const &obj)
+{
+ out << obj.read();
+ return out;
+}
+
+} // namespace fem
diff --git a/cxxfem/src/femTimer.h b/cxxfem/src/femTimer.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c5cf29877adf1268e302735dd732b9583677a1b
--- /dev/null
+++ b/cxxfem/src/femTimer.h
@@ -0,0 +1,32 @@
+#ifndef FEMTIMER_H
+#define FEMTIMER_H
+
+#include "fem.h"
+#include <chrono>
+
+namespace fem
+{
+
+/// a basic timer
+
+class FEM_API Timer
+{
+ std::chrono::high_resolution_clock::time_point t0;
+ double telapsed;
+ bool running;
+
+public:
+ Timer();
+ ~Timer();
+
+ void start();
+ void stop();
+ void reset();
+ double read() const;
+
+ friend FEM_API std::ostream &operator<<(std::ostream &out, Timer const &obj);
+};
+
+} // namespace fem
+
+#endif // FEMTIMER_H
diff --git a/cxxfem/tests/beam2d.py b/cxxfem/tests/beam2d.py
new file mode 100644
index 0000000000000000000000000000000000000000..891407ab21223ee7fbe2f70b43579b0637008b76
--- /dev/null
+++ b/cxxfem/tests/beam2d.py
@@ -0,0 +1,102 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import cxxfem as fem
+import os
+import gmsh
+
+
+def build_square(cx = 0.0, cy = 0.0,
+ Lx = 1.0, Ly = 1.0,
+ nx = 2, ny = 2,
+ structured = True, recombine = True):
+
+ gmsh.model.add("square")
+ lc = 1.0 # no impact
+ gmsh.model.geo.addPoint(cx, cy, 0, lc, 1)
+ gmsh.model.geo.addPoint(cx + Lx, cy, 0, lc, 2)
+ gmsh.model.geo.addPoint(cx + Lx, cy + Ly, 0, lc, 3)
+ gmsh.model.geo.addPoint(cx, cy + Ly, 0, lc, 4)
+
+ gmsh.model.geo.addLine(1, 2, 1)
+ gmsh.model.geo.addLine(2, 3, 2)
+ gmsh.model.geo.addLine(4, 3, 3)
+ gmsh.model.geo.addLine(1, 4, 4)
+
+ gmsh.model.geo.addCurveLoop([1, 2, -3, -4], 1)
+
+ gmsh.model.geo.addPlaneSurface([1], 1)
+
+ gmsh.model.geo.synchronize() # sinon warning "unknown entity..." lors du addPhysicalGroup
+
+ # mesh lines
+ gmsh.model.mesh.setTransfiniteCurve(1, nx + 1)
+ gmsh.model.mesh.setTransfiniteCurve(3, nx + 1)
+ gmsh.model.mesh.setTransfiniteCurve(2, ny + 1)
+ gmsh.model.mesh.setTransfiniteCurve(4, ny + 1)
+
+ # recombine
+ if structured:
+ gmsh.model.mesh.setTransfiniteSurface(1)
+ if recombine:
+ gmsh.model.mesh.setRecombine(2, 1) # quads
+
+ # physical groups - corners
+ pt1 = gmsh.model.addPhysicalGroup(0, [1])
+ pt2 = gmsh.model.addPhysicalGroup(0, [2])
+ pt3 = gmsh.model.addPhysicalGroup(0, [3])
+ pt4 = gmsh.model.addPhysicalGroup(0, [4])
+ gmsh.model.setPhysicalName(0, pt1, "pt1")
+ gmsh.model.setPhysicalName(0, pt2, "pt2")
+ gmsh.model.setPhysicalName(0, pt3, "pt3")
+ gmsh.model.setPhysicalName(0, pt4, "pt4")
+
+ # physical groups - lines
+ bottom = gmsh.model.addPhysicalGroup(1, [1], 1000)
+ right = gmsh.model.addPhysicalGroup(1, [2])
+ top = gmsh.model.addPhysicalGroup(1, [3])
+ left = gmsh.model.addPhysicalGroup(1, [4])
+ gmsh.model.setPhysicalName(1, bottom, "bottom")
+ gmsh.model.setPhysicalName(1, right, "right")
+ gmsh.model.setPhysicalName(1, top, "top")
+ gmsh.model.setPhysicalName(1, left, "left")
+
+ # physical groups - surfaces
+ interior = gmsh.model.addPhysicalGroup(2, [1])
+ gmsh.model.setPhysicalName(2, interior, "interior")
+
+ # finalize and mesh the geometry
+ gmsh.model.mesh.generate(2)
+
+ # gmsh.write("square.msh")
+
+
+if __name__ == "__main__":
+
+ pbl = fem.Problem() #gmsh.initialize()
+ build_square(0.0, 0.0, 10.0, 1.0, 80, 8, True, True)
+ # build_square(0.0, 0.0, 10.0, 1.0, 2, 1, True, True)
+
+ # input()
+
+ m = fem.Medium(pbl, "interior", 100.0, 0.3)
+
+ d1 = fem.Dirichlet(pbl, "left", "X", 0.0)
+ d2 = fem.Dirichlet(pbl, "left", "Y", 0.0)
+ # d3 = fem.Dirichlet(pbl, "pt3", "Y", -0.1) # option 1
+ n3 = fem.Neumann(pbl, "top", "Y", -1e-3) # option 2
+
+ d4 = fem.Dirichlet(pbl, "interior", "Z", 0.0)
+
+
+ solver = fem.Solver(pbl)
+ solver.solve()
+
+ post = fem.Post(solver)
+ post.build_views()
+ post.show_views( [ "stress_tensor", "force_vector" ] )
+ post.deform(5)
+ post.probe("force_vector", 1)
+ post.view()
+
+
\ No newline at end of file
diff --git a/cxxfem/tests/beam3d.py b/cxxfem/tests/beam3d.py
new file mode 100644
index 0000000000000000000000000000000000000000..d6de2a7baf4854357806dfd9ec1804d0e1d1e159
--- /dev/null
+++ b/cxxfem/tests/beam3d.py
@@ -0,0 +1,73 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import cxxfem as fem
+import os
+
+fem.debug_level(0)
+
+if __name__ == "__main__":
+
+ # fine mesh
+ # pars = {
+ # "Lx": 5.0, # width
+ # "Ly": 1.0, # depth
+ # "Lz": 1.0, # height/thickness
+ # "nx": 60,
+ # "ny": 30,
+ # "nz": 30
+ # }
+
+ # same as python
+ pars = {
+ "Lx": 5.0, # width
+ "Ly": 1.0, # depth
+ "Lz": 1.0, # height/thickness
+ "nx": 15,
+ "ny": 2,
+ "nz": 10
+ }
+
+ # mono-element for debug
+ # pars = {
+ # "Lx": 5.0, # width
+ # "Ly": 1.0, # depth
+ # "Lz": 1.0, # height/thickness
+ # "nx": 20,
+ # "ny": 1,
+ # "nz": 1
+ # }
+
+ pbl = fem.Problem()
+ f = os.path.join(os.path.dirname(__file__), 'parallelepiped.geo')
+ pbl.load(f, fem.StringDoubleMap(pars))
+
+ m = fem.Medium(pbl, "Volume", E=100.0, nu=0.3)
+
+ # - mechanical BCs
+
+ # clamping
+ bcs = []
+ for d in ['X','Y','Z']:
+ bcs.append(fem.Dirichlet(pbl, "Left", dof=d, val=0.0))
+
+ # Vertical displacement
+ # bcs.append(fem.Dirichlet(pbl, "Right", dof='Z', val=-1.0)) # prescribed displ
+
+ # bcs.append(fem.Neumann(pbl, "Corner", 'Z', -0.1)) # concentrated force (seems OK)
+ bcs.append(fem.Neumann(pbl, "Top", 'Z', -1e-2)) # distributed load
+ bcs.append(fem.Neumann(pbl, "Top", 'Y', -5e-2)) # distributed load
+
+ solver = fem.Solver(pbl)
+ solver.linear_solver = "pardiso"
+ solver.solve()
+
+ post = fem.Post(solver)
+ # post.build_views(ther=False, smooth=False)
+ post.build_views()
+ post.show_views(['force_vector','stress_tensor'])
+ post.set_camera_3D()
+ # post.write()
+ post.deform()
+ post.view()
+
diff --git a/cxxfem/tests/parallelepiped.geo b/cxxfem/tests/parallelepiped.geo
new file mode 100644
index 0000000000000000000000000000000000000000..e60b38c6da2a4bca6a32d403b5231bee082b483c
--- /dev/null
+++ b/cxxfem/tests/parallelepiped.geo
@@ -0,0 +1,46 @@
+// Parameterized parallelepiped meshed with hexahedra
+
+DefineConstant[ Lx = { 2.0, Min 1, Max 20, Step 1, Name "Lx" } ];// width
+DefineConstant[ Ly = { 2.0, Min 1, Max 20, Step 1, Name "Ly" } ];// depth
+DefineConstant[ Lz = { 2.0, Min 1, Max 20, Step 1, Name "Lz" } ];// height
+DefineConstant[ nx = { 2, Min 1, Max 20, Step 1, Name "nx" } ];// nb of elements along x
+DefineConstant[ ny = { 2, Min 1, Max 20, Step 1, Name "ny" } ];// nb of elements along y
+DefineConstant[ nz = { 2, Min 1, Max 20, Step 1, Name "nz" } ];// nb of elements along z (thickness)
+DefineConstant[
+ recombine = {1, Choices{0="no",1="yes"},
+ Name "recombine"}
+];
+
+d=0.05; // not used
+Point(1) = {0, 0, 0, d};
+Point(2) = {0, Ly, 0, d};
+Point(3) = {Lx, Ly, 0, d};
+Point(4) = {Lx, 0, 0, d};
+Line(1) = {1, 4};
+Line(2) = {4, 3};
+Line(3) = {2, 3};
+Line(4) = {1, 2};
+Line Loop(5) = {4, 3, -2, -1};
+Plane Surface(6) = {5};
+Transfinite Line {4, 2} = ny+1 Using Progression 1;
+Transfinite Line {1, 3} = nx+1 Using Progression 1;
+
+If(recombine == 0)
+ Extrude {0, 0, Lz} { Surface{6}; Layers{nz}; }
+Else
+ Transfinite Surface {6} = {1, 4, 3, 2};
+ Recombine Surface {6};
+ Extrude {0, 0, Lz} { Surface{6}; Layers{nz}; Recombine; }
+EndIf
+
+Physical Volume("Volume") = {1};
+
+Physical Surface("Bottom") = {6}; // Z = 0
+Physical Surface("Top") = {28}; // Z = Lz
+Physical Surface("Front") = {27}; // Y = 0
+Physical Surface("Back") = {19}; // Y = Ly
+Physical Surface("Left") = {15}; // X = 0
+Physical Surface("Right") = {23}; // X = Lx
+
+//+
+Physical Point("Corner", 29) = {4};
diff --git a/envs/cmake.cmd b/envs/cmake.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..ff21a68200b13d9259b7cb708cbf0cc77b180db7
--- /dev/null
+++ b/envs/cmake.cmd
@@ -0,0 +1,3 @@
+@echo off
+::echo cmake.exe -G "MinGW Makefiles" %*
+cmake.exe -G "MinGW Makefiles" %*
\ No newline at end of file
diff --git a/envs/linux-macos.sh b/envs/linux-macos.sh
new file mode 100644
index 0000000000000000000000000000000000000000..9970673d782322100e11d4a020740a14a792be0c
--- /dev/null
+++ b/envs/linux-macos.sh
@@ -0,0 +1,32 @@
+# This file setup an environment which allows you to compile the code
+# on Linux or macOS using the default compiler (gcc or clang)
+#
+# HOW TO USE THIS FILE?
+# open a terminal
+# source ./envs/linux-macos.sh
+# mkdir build
+# cd build
+# cmake ..
+# make
+# [executables are built in the bin/ folder]
+# ctest
+
+# try to get the absolute path of root folder
+DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd && echo x)"; DIR="${DIR%x}"
+abspath=$(dirname $DIR )
+#echo "abspath=${abspath}"
+
+# set the location of gmsh SDK ( **MODIFY THIS LINE FOR YOUR SYSTEM** )
+GMSHSDK=${abspath}/lib/gmsh-sdk
+EIGEN=${abspath}/lib/eigen
+
+# where are gmsh and gmsh-**.so ?
+export PATH=${GMSHSDK}/bin:${GMSHSDK}/lib:${PATH}
+# where is gmsh.h ?
+export INCLUDE=${EIGEN}:${GMSHSDK}/include:${INCLUDE}
+# where is gmsh.lib ?
+export LIB=${GMSHSDK}/lib:${LIB}
+# where is gmsh.py ? (required only if you want to use the python API)
+export PYTHONPATH=${GMSHSDK}/lib:${PYTHONPATH}
+# the following command is only useful for macOS
+export DYLD_LIBRARY_PATH=${GMSHSDK}/lib:${DYLD_LIBRARY_PATH}
diff --git a/envs/make.cmd b/envs/make.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..cd6247a535b99591f68b00b81cc7d587239cb7ab
--- /dev/null
+++ b/envs/make.cmd
@@ -0,0 +1,2 @@
+@echo off
+mingw32-make.exe %*
\ No newline at end of file
diff --git a/envs/win-mingw64.cmd b/envs/win-mingw64.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..b18a4a9309f17f9280ccefc73ab57bee485a3849
--- /dev/null
+++ b/envs/win-mingw64.cmd
@@ -0,0 +1,92 @@
+@echo off
+:: This file opens a terminal which allows you to compile the code with
+:: a 64-bit mingw compiler
+:: https://mingw-w64.org/
+:: https://www.msys2.org/
+::
+:: HOW TO USE THIS FILE?
+:: [check the PATHs below]
+:: [run this file]
+:: mkdir build
+:: cd build
+:: cmake ..
+:: make
+:: ctest
+::
+:: How to clean the "build" folder using cmd line?
+:: cd build
+:: rd /q /s .
+
+echo setting MinGW64 environment...
+
+:: set the location of gmsh SDK / mingw compiler
+set GMSHSDK=%~dp0..\lib\gmsh-sdk
+set EIGEN=%~dp0..\lib\eigen
+set COMPILERPATHS=C:\msys64\mingw64\bin;C:\mingw-w64\mingw64\bin
+
+:: perform some tests...
+IF NOT EXIST %GMSHSDK% (
+ ECHO - gmsh-sdk NOT found in %GMSHSDK%!!
+ PAUSE
+ EXIT /B
+) ELSE (
+ ECHO - gmsh-sdk found in %GMSHSDK%.
+)
+
+IF NOT EXIST %EIGEN% (
+ ECHO - eigen NOT found in %EIGEN%!!
+ PAUSE
+ EXIT /B
+) ELSE (
+ ECHO - eigen found in %EIGEN%.
+)
+
+FOR %%a IN (%COMPILERPATHS%) DO (
+ IF EXIST %%a (
+ ECHO - compiler found in %%a.
+ SET COMPILERPATH=%%a
+ ) ELSE (
+ IF NOT DEFINED COMPILERPATH (
+ ECHO - compiler not found in %%a.
+ )
+ )
+)
+IF NOT DEFINED COMPILERPATH (
+ ECHO =^> compiler NOT found!
+ PAUSE
+ EXIT /B
+)
+
+:: Intel MKL/TBB
+set INTELPATH="C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows"
+IF NOT EXIST %INTELPATH% (
+ ECHO - INTEL libraries NOT found in %INTELPATH%!!
+) ELSE (
+ ECHO - INTEL libraries found.
+ call %INTELPATH%\mkl\bin\mklvars.bat intel64 vs2019
+ call %INTELPATH%\tbb\bin\tbbvars.bat intel64 vs2019
+)
+
+:: where is gmsh.exe and gmsh-**.dll ?
+set PATH=%GMSHSDK%\bin;%GMSHSDK%\lib;%PATH%
+:: where is gmsh.h ?
+set INCLUDE=%EIGEN%;%GMSHSDK%\include;%INCLUDE%
+:: where is gmsh.lib ?
+set LIB=%GMSHSDK%\lib;%LIB%
+:: where is gmsh.py ? (required only if you want to use the python API)
+set PYTHONPATH=%GMSHSDK%\lib;%PYTHONPATH%
+:: add path to the compiler to the PATH
+set PATH=%COMPILERPATH%;%PATH%
+:: add current folder to PATH for cmake/make aliases
+set PATH=%~dp0;%PATH%
+
+:: clear local vars
+set GMSHSDK=
+set EIGEN=
+set COMPILERPATHS=
+set COMPILERPATH=
+
+:: open terminal
+CD /d "%~dp0"
+CD ..
+%comspec% /K
diff --git a/envs/win-msvc.cmd b/envs/win-msvc.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..1aecc53ab4e981839516141e81cb53e434da0026
--- /dev/null
+++ b/envs/win-msvc.cmd
@@ -0,0 +1,68 @@
+@echo off
+:: This file opens a terminal which allows you to compile the code with
+:: Microsoft Visual Studio 2017 (x64) on Windows
+::
+:: HOW TO USE THIS FILE?
+:: [run this file]
+:: mkdir build
+:: cd build
+:: cmake -A x64 ..
+:: cmake --build . --config Release
+:: [executables are built in the bin/Release folder]
+:: ctest -C Release
+
+:: set the location of gmsh SDK / mingw compiler
+set "GMSHSDK=%~dp0..\lib\gmsh-sdk"
+set "EIGEN=%~dp0..\lib\eigen"
+set COMPILERPATH="C:\Program Files (x86)\Microsoft Visual Studio\2019\Community\VC\Auxiliary\Build"
+
+:: perform some tests...
+IF NOT EXIST "%GMSHSDK%" (
+ ECHO - gmsh-sdk NOT found in %GMSHSDK%!!
+ PAUSE
+ EXIT /B
+) ELSE (
+ ECHO - gmsh-sdk found.
+)
+
+IF NOT EXIST "%EIGEN%" (
+ ECHO - eigen NOT found in %EIGEN%!!
+ PAUSE
+ EXIT /B
+) ELSE (
+ ECHO - eigen found.
+)
+
+IF NOT EXIST %COMPILERPATH% (
+ ECHO - compiler NOT found in %COMPILERPATH%!!
+ PAUSE
+ EXIT /B
+) ELSE (
+ ECHO - compiler found.
+)
+
+:: Intel MKL/TBB
+set INTELPATH="C:\Program Files (x86)\IntelSWTools\compilers_and_libraries\windows"
+IF NOT EXIST %INTELPATH% (
+ ECHO - INTEL libraries NOT found in %INTELPATH%!!
+) ELSE (
+ ECHO - INTEL libraries found.
+ call %INTELPATH%\mkl\bin\mklvars.bat intel64 vs2019
+ @REM call %INTELPATH%\tbb\bin\tbbvars.bat intel64 vs2019
+)
+
+:: where is gmsh.exe and gmsh-**.dll ?
+set PATH=%GMSHSDK%\bin;%GMSHSDK%\lib;%PATH%
+:: where is gmsh.h ?
+set INCLUDE=%EIGEN%;%GMSHSDK%\include;%INCLUDE%
+:: where is gmsh.lib ?
+set LIB=%GMSHSDK%\lib;%LIB%
+:: where is gmsh.py ? (required only if you want to use the python API)
+set PYTHONPATH=%GMSHSDK%\lib;%PYTHONPATH%
+
+:: set the environment of the msvc compiler
+CD /d "%~dp0"
+CD ..
+:: problem with mkl_rt (mkl_intel_thread cannot be loaded)
+:: %comspec% /k "%COMPILERPATH%\vcvars64.bat"
+%comspec% /k
diff --git a/lib/.gitignore b/lib/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..67a64465eb935decc500dbb91de7b25924e2fbca
--- /dev/null
+++ b/lib/.gitignore
@@ -0,0 +1,2 @@
+eigen/
+gmsh-sdk/
\ No newline at end of file
diff --git a/lib/README.md b/lib/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..ed65683f8679789ff6fd7723c78a9217eb4bdf03
--- /dev/null
+++ b/lib/README.md
@@ -0,0 +1,3 @@
+# ./lib
+
+This folder contains scripts for downloading the 2 external libraries that are required for this project.
\ No newline at end of file
diff --git a/lib/get_eigen.cmd b/lib/get_eigen.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..cf848c53f6846dd7f9f005c456875ef7568d1ef0
--- /dev/null
+++ b/lib/get_eigen.cmd
@@ -0,0 +1,18 @@
+::@echo off
+:: download eigen & extract to current folder
+
+set version=3.3.9
+set file=eigen-%version%.zip
+
+:: download file
+del %file% >nul 2>&1
+bitsadmin /transfer get_eigen /dynamic /download /priority foreground "https://gitlab.com/libeigen/eigen/-/archive/%version%/%file%" "%CD%\%file%"
+
+:: unzip
+rd /Q/S eigen-%version% >nul 2>&1
+powershell.exe -nologo -noprofile -command "& { Add-Type -A 'System.IO.Compression.FileSystem'; [IO.Compression.ZipFile]::ExtractToDirectory('%CD%\%file%', '%CD%'); }"
+del %file% >nul 2>&1
+
+:: rename folder
+rd /Q/S eigen >nul 2>&1
+ren eigen-%version% eigen
diff --git a/lib/get_eigen.sh b/lib/get_eigen.sh
new file mode 100644
index 0000000000000000000000000000000000000000..aea96329790dd7df0f5628c80ec231f3ea72b071
--- /dev/null
+++ b/lib/get_eigen.sh
@@ -0,0 +1,11 @@
+#!/bin/bash
+
+set -e
+
+VERSION=3.3.9
+
+wget -q https://gitlab.com/libeigen/eigen/-/archive/${VERSION}/eigen-${VERSION}.tar.gz
+tar xzf eigen-${VERSION}.tar.gz
+rm eigen-${VERSION}.tar.gz
+rm -rf eigen
+mv eigen-${VERSION} eigen
diff --git a/lib/get_gmsh.cmd b/lib/get_gmsh.cmd
new file mode 100644
index 0000000000000000000000000000000000000000..d9a71d6be81d90b6311bfdec9a43b89896f25931
--- /dev/null
+++ b/lib/get_gmsh.cmd
@@ -0,0 +1,28 @@
+::@echo off
+:: download gmsh-sdk & extract to current folder
+
+set version=4.10.2
+
+set gmsh=gmsh-%version%-Windows64-sdk
+set file=%gmsh%.zip
+
+:: download file
+del %file% >nul 2>&1
+bitsadmin /transfer get_gmsh /dynamic /download /priority foreground "https://gmsh.info/bin/Windows/%file%" "%CD%\%file%"
+
+:: unzip
+rd /Q/S %gmsh% >nul 2>&1
+powershell.exe -nologo -noprofile -command "& { Add-Type -A 'System.IO.Compression.FileSystem'; [IO.Compression.ZipFile]::ExtractToDirectory('%CD%\%file%', '%CD%'); }"
+del %file% >nul 2>&1
+
+:: rename folder
+rd /Q/S gmsh-sdk >nul 2>&1
+ren %gmsh% gmsh-sdk
+
+:: copy dll to bin folder so that double-clic on gmsh.exe works!
+move /Y gmsh-sdk\lib\gmsh-*.dll gmsh-sdk\bin\
+:: IMPORTANT NOTE: do not leave gmsh-*.dll next to *.py because gmsh.py
+:: will try to load this one first.
+:: However any C++ program will prefer gmsh-*.dll in the bin folder
+:: => any combined python/C++ program will load both copies in memory
+:: => C++ calls and python call will not operate on the same memory.
\ No newline at end of file
diff --git a/lib/get_gmsh.sh b/lib/get_gmsh.sh
new file mode 100644
index 0000000000000000000000000000000000000000..4b8cb825426bbfbfe05520be1e48f8f1f058ab0c
--- /dev/null
+++ b/lib/get_gmsh.sh
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+set -e
+
+# VERSION=4.9.3
+VERSION=4.10.2
+
+if [[ "$OSTYPE" == "linux-gnu" ]]; then
+ wget -q https://gmsh.info/bin/Linux/gmsh-${VERSION}-Linux64-sdk.tgz
+ tar xzf gmsh-${VERSION}-Linux64-sdk.tgz
+ rm gmsh-${VERSION}-Linux64-sdk.tgz
+ rm -rf gmsh-sdk
+ mv gmsh-${VERSION}-Linux64-sdk gmsh-sdk
+elif [[ "$OSTYPE" == "darwin"* ]]; then
+ wget -q https://gmsh.info/bin/MacOSX/gmsh-${VERSION}-MacOSX-sdk.tgz
+ tar xzf gmsh-${VERSION}-MacOSX-sdk.tgz
+ rm gmsh-${VERSION}-MacOSX-sdk.tgz
+ rm -rf gmsh-sdk
+ mv gmsh-${VERSION}-MacOSX-sdk gmsh-sdk
+else
+ echo "unknown system"
+fi
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
new file mode 100644
index 0000000000000000000000000000000000000000..6f8e08fd56c0a8af11b4c43c84e3fdcf8998ac3c
--- /dev/null
+++ b/models/bonemodel2.py
@@ -0,0 +1,324 @@
+# -*- coding: utf-8 -*-
+# Bone FEA model (formely: "mandible model")
+#
+# main cxxfem script with default parameters.
+# R.Boman
+
+import cxxfem as fem
+from . import boneload
+import os
+import gmsh
+import numpy as np
+
+
+def parms(d={}):
+ """default parameters of the model
+ """
+ p = {}
+
+ # Geomagic/MeshLab input files
+ path = 'dolicorhynchops/10k'
+ # mandible mesh (.ply/.geo/.msh)
+ p['bone'] = f'{path}/mandible.ply'
+ # one or several vertices (.ply/.off) (used if p['food']=='fixteeth')
+ p['contact_pts'] = f'{path}/teeth.off'
+ # left temporomandibular joint - 1 vertex (.ply/.off/coordinates)
+ p['axis_pt1'] = [0.1458893, -73.13895, 227.3909]
+ # right temporomandibular joint - 1 vertex (.ply/.off/coordinates)
+ p['axis_pt2'] = f'{path}/RTMJ.off'
+
+ p['muscles'] = [
+ {
+ 'file': f'{path}/Lmuscle.ply',
+ 'force': 100.,
+ 'focalpt': f'{path}/LmuscleF.off',
+ 'method': 'T' # loading model: 'U', 'T', 'T+N'
+ },
+ {
+ 'file': f'{path}/Rmuscle.off',
+ 'force': 100.,
+ 'focalpt': f'{path}/RmuscleF.off',
+ 'method': 'T'
+ }
+ ]
+ p['food'] = 'fixteeth' # type of food: 'fixteeth' or 'cylinder'
+ p['fixations'] = {
+ 'contact_pts': ['x', 'y', 'z'], # constrained degrees of freedom
+ 'axis_pt1': ['x', 'y', 'z'],
+ 'axis_pt2': ['x', 'y', 'z']
+ }
+ # rigid cylinder (used if p['food']=='cylinder')
+ R = 20
+ p['cylinder'] = {
+ 'origin': [-R-8, 0, -150],
+ 'radius': R,
+ 'width': 100,
+ 'friction': 0.1,
+ 'penalty': 1e3
+ }
+
+ # material properties
+ p['density'] = 1.850e-9 # [T/mm³] - bone: 1.850 kg/l
+ p['Young'] = 17000. # [MPa] elastic modulus - bone: 17-20 GPa
+ p['Poisson'] = 0.3 # [-] Poisson's ratio
+
+ # numerical parameters
+ p['tolNR'] = 1e-6 # [-] equilibrium tolerance
+ p['dt0'] = 1.0 # [s] time step size
+
+ # gmsh toolbox
+ p['use_gmshOld'] = False # use old gmsh interface
+
+ p.update(d)
+ return p
+
+
+def solve(p={}):
+
+ p = parms(p)
+
+ import json
+ # print(f'parameters={json.dumps(p, indent=4, sort_keys=True)}')
+
+ pbl = fem.Problem()
+
+ gmsh.initialize()
+ gmsh.option.setNumber("General.Terminal", 1)
+ # load main mandible mesh (volume or surface - it will be meshed in 3D)
+ import_mesh(p['bone'])
+
+ # define groups
+
+ # group (3,1)='all' is the group containing all volume elements
+ groups = {}
+ groups['all'] = (3, 1)
+ nodeTags, coord = gmsh.model.mesh.getNodesForPhysicalGroup(*groups['all'])
+ print(f"the mandible mesh has {len(nodeTags)} nodes")
+
+ # extract external surface of the mesh => groups['surf']
+ # # nods_no: tags of the nodes
+ # # nods_pos: coordinates of these nodes
+ groups['surf'] = (2, 1)
+ nods_no, nods_pos = gmsh.model.mesh.getNodesForPhysicalGroup(
+ *groups['surf'])
+ print(f"the mandible surface has {len(nods_no)} nodes")
+ nods_pos = np.reshape(nods_pos, (len(nods_no), 3))
+ # print(f'nods_pos={nods_pos}')
+
+ # load other surface meshes => groups['axis_pt1'], ...
+ for meshname in ['axis_pt1', 'axis_pt2']:
+ create_group(p[meshname], nods_no, nods_pos, groups, meshname)
+
+ # load muscle groups
+ mgroups = {} # stores muscle group data and loads
+ for muscle in p['muscles']:
+ # load focal point
+ if isinstance(muscle['focalpt'], str):
+ fullpath = os.path.join(
+ os.path.dirname(__file__), muscle['focalpt'])
+ focalnodes, _ = boneload.load_msh(fullpath)
+ else: # coordinates in array or tuple
+ focalnodes = [muscle['focalpt']]
+
+ # load surface mesh => groups[name (from filename)]
+ name, nodes, tris, ntags = \
+ create_group(muscle['file'], nods_no, nods_pos, groups)
+ # create loads on the surface
+ loads = boneload.create_loads(nodes, tris,
+ total_force=muscle['force'], target=focalnodes[0], method=muscle['method'])
+ # store data in a structure
+ mgroups[name] = MuscleGroup(nodes, tris, ntags, loads)
+
+ if p['food'] == 'fixteeth':
+ # teeth are fixed along chosen directions
+ create_group(p['contact_pts'], nods_no,
+ nods_pos, groups, 'contact_pts')
+
+
+ # --------------------------------------------------------------------------
+ print('== ADDING MEDIUM...')
+
+ do_not_delete = []
+
+ # material
+ do_not_delete.append( fem.Medium(pbl, "all", E=p['Young'], nu=p['Poisson']) )
+
+ print('== ADDING DIRICHLETs...')
+ # axis of rotation
+ for gname in ['axis_pt1', 'axis_pt2']:
+ for d in p['fixations'][gname]:
+ do_not_delete.append(fem.Dirichlet(pbl, gname, dof=d.upper(), val=0.0) )
+
+ print('== ADDING NEUMANNs...')
+ # mshpts = geometry.getMesh().getPointSet()
+ for name, gr in mgroups.items():
+ # print(f'len(gr.ntags)={len(gr.ntags)}')
+ for i,load in enumerate(gr.loads):
+ do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'X', load.x[0]) ) # TODO: regler pbl de cast "gr.ntags[i] => size_t"
+ do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'Y', load.x[1]) ) # TODO: regler "dof=" qui marche pas
+ do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'Z', load.x[2]) )
+
+ # type of food
+
+ if p['food'] == 'fixteeth':
+ # teeth are fixed along chosen directions
+ for d in p['fixations']['contact_pts']:
+ do_not_delete.append(fem.Dirichlet(pbl, 'contact_pts', dof=d.upper(), val=0.0) )
+ else:
+ raise Exception('food = "fixteeth" only supported')
+
+
+ # python only!
+ # do_not_delete.append(fem.Dirichlet(pbl, "all", dof='T', val=0.0) )
+
+
+ # Time integration scheme
+ # print (pbl)
+ print('SOLVE.......')
+
+ solver = fem.Solver(pbl)
+ solver.linear_solver = "pardiso"
+ solver.solve()
+
+ post = fem.Post(solver)
+ post.build_views()
+ post.show_views([ 'smooth_stress_tensor' ])
+ post.set_camera_3D()
+ post.write()
+ post.deform()
+
+ print('extracting results...')
+ print(f'\tinternal energy = {solver.int_energy:.3f} N.mm')
+ for grpname in ['axis_pt1', 'axis_pt2', 'contact_pts']:
+ v = np.array(post.probe('force_vector', grpname))
+ v = np.reshape(v, (v.size//3, 3))
+ v = np.sum(v, axis=0)
+ print(f'\t{grpname}_Fx = {v[0]:.3f} N')
+ print(f'\t{grpname}_Fy = {v[1]:.3f} N')
+ print(f'\t{grpname}_Fz = {v[2]:.3f} N')
+
+ for name, group in mgroups.items():
+ v=np.array( [0.0, 0.0, 0.0] )
+ for tag in group.ntags:
+ v += np.array(post.probe('force_vector', int(tag)))
+ print(f'\t{name}_Fx = {v[0]:.3f} N')
+ print(f'\t{name}_Fy = {v[1]:.3f} N')
+ print(f'\t{name}_Fz = {v[2]:.3f} N')
+
+ post.view()
+
+
+# ------------------------------------------------------------------------------
+# utility functions involing Gmsh objects
+
+
+def import_mesh(filename):
+ """import the volume mesh
+ input file type can be:
+ - .msh: the mesh is simply loaded
+ - .geo: the mesh is generated
+ - .ply/.stl: a .geo is written in the workspace and loaded into gmsh
+ """
+ # import gmsh file (accept .geo, .msh, .ply or .stl)
+ mandible_gmsh = filename
+ # make the path absolute
+ if not os.path.isabs(mandible_gmsh):
+ mandible_gmsh = os.path.join(os.path.dirname(__file__), mandible_gmsh)
+ # check whether file exists
+ if not os.path.isfile(mandible_gmsh):
+ raise Exception(f'{mandible_gmsh} not found!')
+ # extract extension
+ noext, ext = os.path.splitext(mandible_gmsh)
+ if ext in ['.ply', '.stl']:
+ plyfile = mandible_gmsh
+ print(f'loading file {plyfile} into Gmsh...')
+ # gmsh.model.add("fossil")
+ gmsh.merge(plyfile)
+ gmsh.model.geo.addSurfaceLoop(surfaceTags=[1], tag=1)
+ gmsh.model.geo.addVolume(shellTags=[1], tag=1)
+ gmsh.model.geo.addPhysicalGroup(dim=2, tags=[1], tag=1)
+ gmsh.model.setPhysicalName(dim=2, tag=1, name='surf')
+ gmsh.model.geo.addPhysicalGroup(dim=3, tags=[1], tag=1)
+ gmsh.model.setPhysicalName(dim=3, tag=1, name='all')
+ elif ext in ['.geo', '.msh']:
+ gmsh.merge(mandible_gmsh)
+ else:
+ raise Exception(
+ f'Unknown extension: {ext}, please use .ply, .stl, .msh or .geo')
+
+ # generate mesh
+ gmsh.model.geo.synchronize()
+ gmsh.model.mesh.generate(3)
+ # gmsh.fltk.run()
+
+
+def create_group(mesh_file_or_coord, all_no, all_coords, groups, gname=None):
+ """loads the mesh file "mesh_file_or_coord" and identify vertices among vertices
+ stored in all_coords.
+ Then, the routine creates a group named "gname" (built from mesh_file_or_coord if None)
+
+ mesh_file_or_coord can be either:
+ a filename
+ an array of points: [ [x1,y1,z1], [x2,y2,z2], ... ]
+ one single point: [x1,y1,z1]
+ """
+ if isinstance(mesh_file_or_coord, str):
+ # load mesh file
+ print(f'\n* create_group {mesh_file_or_coord}...')
+ fullpath = os.path.join(os.path.dirname(__file__), mesh_file_or_coord)
+ nodes, tris = boneload.load_msh(fullpath)
+ elif type(mesh_file_or_coord) in [list, tuple]:
+ print(f'\n* create_group {gname}...')
+ if type(mesh_file_or_coord[0]) in [list, tuple]:
+ # array of points
+ nodes = mesh_file_or_coord
+ else:
+ # one single point
+ nodes = [mesh_file_or_coord]
+ print(f'nodes={nodes}')
+ print(f'mesh_file_or_coord={mesh_file_or_coord}')
+ tris = []
+ else:
+ raise Exception(f'Data is not understood: {mesh_file_or_coord}')
+ # identify nodes of the mesh among the given vertices of the main mesh
+ # ntags = boneload.identify_nodes_brutal(nodes, all_no, all_coords)
+ ntags = boneload.identify_nodes_SaP(nodes, all_no, all_coords)
+ # print(f'ntags={ntags}')
+
+ if gname is None:
+ gname, _ = os.path.splitext(os.path.basename(mesh_file_or_coord))
+
+ # create a Gmsh group
+ # => a discrete entity with "points" elements
+ print(f'** adding gmsh group "{gname}"...')
+ etag = gmsh.model.addDiscreteEntity(dim=0)
+ point_type = gmsh.model.mesh.getElementType("Point", 1)
+ # print(f'point_type={point_type}')
+ gmsh.model.mesh.addElementsByType(etag, point_type, [], ntags)
+ ptag = gmsh.model.geo.addPhysicalGroup(dim=0, tags=[etag])
+ gmsh.model.setPhysicalName(dim=0, tag=ptag, name=gname)
+
+ groups[gname] = (0, ptag)
+
+ # required for the nodes to be linked to the physical group
+ gmsh.model.geo.synchronize()
+ nods, pos = gmsh.model.mesh.getNodesForPhysicalGroup(*groups[gname])
+ # print(f'nods={nods}')
+ # print(f'pos={pos}')
+
+ print(f'group {groups[gname]} ({gname}) '
+ f'created from {mesh_file_or_coord} ({len(nods)} nodes)')
+
+ return gname, nodes, tris, ntags
+
+
+class MuscleGroup:
+ """convenient structure for storing muscle data and loads
+ """
+
+ def __init__(self, nodes, tris, ntags, loads):
+ self.nodes = nodes
+ self.tris = tris
+ self.ntags = ntags
+ self.loads = loads
+
diff --git a/models/dolicorhynchops/dolicorhynchops_10k.py b/models/dolicorhynchops/dolicorhynchops_10k.py
index 6087c273a6ee9a90004106a5cb195006791338af..0ebd8da4ff458f3eb5fe53ccc87792cdf50011d5 100644
--- a/models/dolicorhynchops/dolicorhynchops_10k.py
+++ b/models/dolicorhynchops/dolicorhynchops_10k.py
@@ -9,9 +9,9 @@ def parms(d={}):
path = 'dolicorhynchops/10k'
p['bone'] = f'{path}/mandible.stl'
p['contact_pts'] = [[-10.20362, -17.46838, -229.9061],
- [-11.92466, 26.3042, -229.5354]]
+ [-11.92466, 26.3042, -229.5354]]
p['axis_pt1'] = f'{path}/LTMJ.off'
- p['axis_pt2'] = [ -8.716309, 79.13171, 233.8385 ]
+ p['axis_pt2'] = [-8.716309, 79.13171, 233.8385]
p['muscles'] = [
{
'file': f'{path}/Lmuscle.stl',
@@ -45,3 +45,8 @@ def parms(d={}):
def getMetafor(p={}):
import bonemodel as model
return model.getMetafor(parms(p))
+
+
+if __name__ == "__main__":
+ import models.bonemodel2 as model
+ model.solve(parms())
diff --git a/run.py b/run.py
new file mode 100644
index 0000000000000000000000000000000000000000..4bb930750261ec4f32ba0b32c8c1a60cd5f0da4d
--- /dev/null
+++ b/run.py
@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# runs a test as if it was installed
+
+import os, platform
+
+if 'Windows' in platform.uname():
+ # remove tbb from the PATH
+ # otherwise: "Intel MKL FATAL ERROR: Cannot load mkl_intel_thread.dll"
+ paths = os.environ['PATH'].split(';')
+ corr_path = [ p for p in paths if not 'tbb' in p]
+ os.environ['PATH'] = ';'.join(corr_path)
+ for p in corr_path:
+ print(p)
+
+if __name__ == "__main__":
+ import sys
+ import os
+ # adds "." to the pythonpath
+ thisdir = os.path.split(os.path.abspath(__file__))[0]
+ thisdir = os.path.normcase(thisdir)
+ print(f'adding {thisdir} to pythonpath')
+ sys.path.append(thisdir)
+
+ # add binary dir to PYTHONPATH
+ pyexe = os.path.basename(sys.executable)
+ print(f'pyexe={pyexe}')
+ sys.path.append(os.path.join(thisdir, 'cxxfem',
+ 'build', 'bin')) # gcc/mingw
+ sys.path.append(os.path.join(thisdir, 'cxxfem',
+ 'build', 'bin', 'Release')) # msvc
+
+ # display env variables
+ try:
+ print(f"OMP_NUM_THREADS={os.environ['OMP_NUM_THREADS']}")
+ except:
+ print(f"OMP_NUM_THREADS=[not set]")
+
+ # parse args
+ import argparse
+ parser = argparse.ArgumentParser()
+ parser.add_argument(
+ "-v", "--verb", help="increase output verbosity", action="count", default=0)
+ parser.add_argument("--nogui", help="disable any graphical output",
+ action="store_true")
+ parser.add_argument('file', nargs='*', help='python files')
+ args = parser.parse_args()
+
+ # run all tests sequentially
+ for testname in args.file:
+ testname = os.path.abspath(testname)
+ testname = os.path.normcase(testname) # F:/ => f:/ on Windows
+ # create workspace
+ common = os.path.commonprefix((testname, thisdir + os.sep))
+ resdir = testname[len(common):].replace(os.sep, "_")
+ resdir, ext = os.path.splitext(resdir)
+ wdir = os.path.join('workspace', resdir)
+ print('workspace=', wdir)
+ if not os.path.isdir(wdir):
+ os.makedirs(wdir)
+ os.chdir(wdir)
+ if ext == '.py':
+ # python script
+ __file__ = testname
+ exec(open(testname, encoding='utf-8').read())
diff --git a/run.spec b/run.spec
new file mode 100644
index 0000000000000000000000000000000000000000..d795a74140de6545cc5f3d3ba26e24df71112b4d
--- /dev/null
+++ b/run.spec
@@ -0,0 +1,60 @@
+# -*- mode: python ; coding: utf-8 -*-
+
+# pyinstaller --noconfirm run.spec
+
+block_cipher = None
+
+
+a = Analysis(
+ ['run.py'],
+ pathex=['cxxfem\\build\\bin\\Release'],
+ binaries=[
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_vml_avx512.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_rt.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_intel_thread.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_core.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_avx512.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_avx.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\mkl\mkl_def.dll', '.'),
+ (r'C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2019.5.281\windows\redist\intel64_win\compiler\libiomp5md.dll', '.')
+ ],
+ datas=[],
+ hiddenimports=['cxxfem','femi','gmsh','bonemodel'],
+ hookspath=[],
+ hooksconfig={},
+ runtime_hooks=[],
+ excludes=[],
+ win_no_prefer_redirects=False,
+ win_private_assemblies=False,
+ cipher=block_cipher,
+ noarchive=False,
+)
+pyz = PYZ(a.pure, a.zipped_data, cipher=block_cipher)
+
+exe = EXE(
+ pyz,
+ a.scripts,
+ [],
+ exclude_binaries=True,
+ name='run',
+ debug=False,
+ bootloader_ignore_signals=False,
+ strip=False,
+ upx=True,
+ console=True,
+ disable_windowed_traceback=False,
+ argv_emulation=False,
+ target_arch=None,
+ codesign_identity=None,
+ entitlements_file=None,
+)
+coll = COLLECT(
+ exe,
+ a.binaries,
+ a.zipfiles,
+ a.datas,
+ strip=False,
+ upx=True,
+ upx_exclude=[],
+ name='run',
+)