add C++ FEM solver

.clang-format

+---
+# 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
+...

.editorconfig

+# EditorConfig is awesome: https://EditorConfig.org
+
+root = true
+
+[*]
+charset = utf-8
+indent_style = space
+indent_size = 4

.gitignore

@@ -8,3 +8,4 @@ __pycache__/
 *.zip
 *.7z
 *.rar
+build/

cxxfem/CMakeLists.txt

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

cxxfem/__init__.py

+# -*- coding: utf-8 -*-
+# fem MODULE initialization file
+
+from femi import *

cxxfem/_src/CMakeLists.txt

+# 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})

cxxfem/_src/femi.i

+// 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"

cxxfem/config.h.in

+#cmakedefine PROJECT_SOURCE_DIR "@PROJECT_SOURCE_DIR@"
+#cmakedefine FEM_USE_MKL

cxxfem/exe/CMakeLists.txt

+
+# 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()

cxxfem/exe/beam.cpp

+// 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;
+    }
+}

cxxfem/exe/build_square.h

+#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

cxxfem/exe/simple_shear.cpp

+// 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;
+    }
+}

cxxfem/exe/square.cpp

+// 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;
+    }
+}

cxxfem/src/CMakeLists.txt

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

cxxfem/src/OpenMP.h

+#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

cxxfem/src/fem.cpp

+#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

cxxfem/src/fem.h

+#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

cxxfem/src/femDOF.cpp

+#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

cxxfem/src/femDOF.h

+#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

cxxfem/src/femDirichlet.cpp

+#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