diff --git a/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake b/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake
index eafef3917..aa3e071be 100644
--- a/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake
+++ b/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake
@@ -1,820 +1,822 @@
 #===============================================================================
 # @file   CMakePackagesSystem.cmake
 #
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @brief  Set of macros used by the package system, internal functions
 #
 # @section LICENSE
 #
 # Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # Akantu is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Akantu. If not, see <http://www.gnu.org/licenses/>.
 #
 #===============================================================================
 
 # ==============================================================================
 # "Private" Accessors
 # ==============================================================================
 # ------------------------------------------------------------------------------
 # Real name
 # ------------------------------------------------------------------------------
 function(_package_get_real_name pkg_name real_name)
   set(${real_name} ${${pkg_name}} PARENT_SCOPE)
 endfunction()
 
 function(_package_set_real_name pkg_name real_name)
   set(${pkg_name} ${real_name} CACHE INTERNAL "" FORCE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Option name
 # ------------------------------------------------------------------------------
 function(_package_declare_option pkg_name)
   string(TOUPPER "${PROJECT_NAME}" _project)
   _package_get_real_name(${pkg_name} _real_name)
   string(TOUPPER "${_real_name}" _u_package)
 
   _package_get_nature(${pkg_name} _nature)
 
   if(${_nature} MATCHES "internal" OR ${_nature} MATCHES "meta")
     set(_opt_name ${_project}_${_u_package})
   elseif(${_nature} MATCHES "external")
     set(_opt_name ${_project}_USE_${_u_package})
   else()
     set(_opt_name UNKNOWN_NATURE_${_project}_${_u_package})
   endif()
 
   _package_set_variable(OPTION_NAME ${pkg_name} ${_opt_name})
 endfunction()
 
 function(_package_get_option_name pkg_name opt_name)
   _package_get_variable(OPTION_NAME ${pkg_name} _opt_name)
   set(${opt_name} ${_opt_name} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Set if system package or compile external lib
 # ------------------------------------------------------------------------------
 function(_package_set_system_option pkg_name default)
   string(TOUPPER "${PROJECT_NAME}" _project)
   _package_get_real_name(${pkg_name} _real_name)
   string(TOUPPER "${_real_name}" _u_package)
 
   option(${_project}_USE_SYSTEM_${_u_package}
     "Should akantu compile the third-party: \"${_real_name}\"" ${default})
   mark_as_advanced(${_project}_USE_SYSTEM_${_u_package})
 endfunction()
 
 function(_package_use_system pkg_name use)
   string(TOUPPER "${PROJECT_NAME}" _project)
   _package_get_real_name(${pkg_name} _real_name)
   string(TOUPPER "${_real_name}" _u_package)
   if(DEFINED ${_project}_USE_SYSTEM_${_u_package})
     set(${use} ${${_project}_USE_SYSTEM_${_u_package}} PARENT_SCOPE)
   else()
     set(${use} TRUE PARENT_SCOPE)
   endif()
 endfunction()
 
 function(_package_set_system_script pkg_name script)
   _package_set_variable(COMPILE_SCRIPT ${pkg_name} "${script}")
 endfunction()
 
 function(_package_add_third_party_script_variable pkg_name var)
   _package_set_variable(VARIABLE_${var} ${pkg_name} "${ARGN}")
   set(${var} ${ARGN} PARENT_SCOPE)
 endfunction()
 
 function(_package_load_third_party_script pkg_name)
   if(${pkg_name}_COMPILE_SCRIPT)
     # set the stored variable
     get_cmake_property(_all_vars VARIABLES)
     foreach(_var ${_all_vars})
       if(_var MATCHES "^${pkg_name}_VARIABLE_.*")
 	string(REPLACE "${pkg_name}_VARIABLE_" "" _orig_var "${_var}")
 	set(${_orig_var} ${${_var}})
       endif()
     endforeach()
 
     _package_get_real_name(${pkg_name} _name)
     string(TOUPPER "${_name}" _u_name)
 
     _package_get_option_name(${pkg_name} _opt_name)
     if(${_opt_name}_VERSION)
       set(_version " (version ${${_opt_name}_VERSION})")
     elseif(${_u_name}_VERSION)
       set(_version " (version ${${_u_name}_VERSION})")
     endif()
 
     # load the script
     message(STATUS "${_name}: building as third-party${_version}")
     include(ExternalProject)
     include(${${pkg_name}_COMPILE_SCRIPT})
   endif()
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Nature
 # ------------------------------------------------------------------------------
 function(_package_set_nature pkg_name nature)
   _package_set_variable(NATURE ${pkg_name} ${nature})
 endfunction()
 
 function(_package_get_nature pkg_name nature)
   _package_get_variable(NATURE ${pkg_name} _nature "unknown")
   set(${nature} ${_nature} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Description
 # ------------------------------------------------------------------------------
 function(_package_set_description pkg_name desc)
   _package_set_variable(DESC ${pkg_name} ${desc})
 endfunction()
 
 function(_package_get_description pkg_name desc)
   _package_get_variable(DESC ${pkg_name} _desc "No description set for the package ${${pkg_name}} (${pkg_name})")
   set(${desc} ${_desc} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Package file name
 # ------------------------------------------------------------------------------
 function(_package_set_filename pkg_name file)
   _package_set_variable(FILE ${pkg_name} ${file})
 endfunction()
 
 function(_package_get_filename pkg_name file)
   _package_get_variable(FILE ${pkg_name} _file "No filename set for the package ${${pkg_name}}")
   set(${file} ${_file} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Source folder
 # ------------------------------------------------------------------------------
 function(_package_set_sources_folder pkg_name src_folder)
   _package_set_variable(SRC_FOLDER ${pkg_name} ${src_folder})
 endfunction()
 
 function(_package_get_sources_folder pkg_name src_folder)
   _package_get_variable(SRC_FOLDER ${pkg_name} _src_folder)
   set(${src_folder} ${_src_folder} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Test folder
 # ------------------------------------------------------------------------------
 function(_package_set_tests_folder pkg_name test_folder)
   _package_set_variable(TEST_FOLDER ${pkg_name} ${test_folder})
 endfunction()
 
 function(_package_get_tests_folder pkg_name test_folder)
   _package_get_variable(TEST_FOLDER ${pkg_name} _test_folder)
   set(${test_folder} ${_test_folder} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Manual folder
 # ------------------------------------------------------------------------------
 function(_package_set_manual_folder pkg_name manual_folder)
   _package_set_variable(MANUAL_FOLDER ${pkg_name} ${manual_folder})
 endfunction()
 
 function(_package_get_manual_folder pkg_name manual_folder)
   _package_get_variable(MANUAL_FOLDER ${pkg_name} _manual_folder)
   set(${manual_folder} ${_manual_folder} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Extra option for the find_package
 # ------------------------------------------------------------------------------
 function(_package_set_find_package_extra_options pkg_name)
   _package_set_variable(FIND_PKG_OPTIONS ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_find_package_extra_options pkg_name options)
   _package_get_variable(FIND_PKG_OPTIONS ${pkg_name} _options)
   set(${options} ${_options} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Compilation flags
 # ------------------------------------------------------------------------------
 function(_package_set_compile_flags pkg_name)
   _package_set_variable(COMPILE_FLAGS ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_compile_flags pkg_name flags)
   _package_get_variable(COMPILE_FLAGS ${pkg_name} _flags)
   set(${flags} ${_flags} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Include dir
 # ------------------------------------------------------------------------------
 function(_package_set_include_dir pkg_name)
   _package_set_variable(INCLUDE_DIR ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_include_dir pkg_name include_dir)
   _package_get_variable(INCLUDE_DIR ${pkg_name} _include_dir "")
   set(${include_dir} ${_include_dir} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Libraries
 # ------------------------------------------------------------------------------
 function(_package_set_libraries pkg_name)
   _package_set_variable(LIBRARIES ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_libraries pkg_name libraries)
   _package_get_variable(LIBRARIES ${pkg_name} _libraries "")
   set(${libraries} ${_libraries} PARENT_SCOPE)
 endfunction()
 
 function(_package_add_libraries pkg_name)
   _package_add_to_variable(LIBRARIES ${pkg_name} ${ARGN})
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Extra dependencies like custom commands of ExternalProject
 # ------------------------------------------------------------------------------
 function(_package_add_extra_dependency pkg_name)
   _package_add_to_variable(EXTRA_DEPENDENCY ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_rm_extra_dependency pkg_name dep)
   _package_remove_from_variable(EXTRA_DEPENDENCY ${pkg_name} ${dep})
 endfunction()
 
 function(_package_set_extra_dependencies pkg)
   _package_set_variable(EXTRA_DEPENDENCY ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_extra_dependencies pkg deps)
   _package_get_variable(EXTRA_DEPENDENCY ${pkg_name} _deps "")
   set(${deps} ${_deps} PARENT_SCOPE)
 endfunction()
 
 function(_package_unset_extra_dependencies pkg_name)
   _package_variable_unset(EXTRA_DEPENDENCY ${pkg_name})
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Activate/deactivate
 # ------------------------------------------------------------------------------
 function(_package_activate pkg_name)
   _package_set_variable(STATE ${pkg_name} ON)
 endfunction()
 
 function(_package_deactivate pkg_name)
   _package_set_variable(STATE ${pkg_name} OFF)
 endfunction()
 
 function(_package_is_activated pkg_name act)
   _package_get_variable(STATE ${pkg_name} _state OFF)
   if(_state)
     set(${act} TRUE PARENT_SCOPE)
   else()
     set(${act} FALSE PARENT_SCOPE)
   endif()
 endfunction()
 
 function(_package_is_deactivated pkg_name act)
   _package_get_variable(STATE ${pkg_name} _state OFF)
   if(NOT _state)
     set(${act} TRUE PARENT_SCOPE)
   else()
     set(${act} FALSE PARENT_SCOPE)
   endif()
 endfunction()
 
 function(_package_unset_activated pkg_name)
   _package_variable_unset(STATE ${pkg_name})
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Export list
 # ------------------------------------------------------------------------------
 function(_package_add_to_export_list pkg_name)
   _package_add_to_variable(EXPORT_LIST ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_export_list pkg_name export_list)
   _package_get_variable(EXPORT_LIST ${pkg_name} _export_list)
   set(${export_list} ${_export_list} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Direct dependencies
 # ------------------------------------------------------------------------------
 function(_package_add_dependencies pkg_name)
   _package_add_to_variable(DEPENDENCIES ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_dependencies pkg_name dependencies)
   _package_get_variable(DEPENDENCIES ${pkg_name} _dependencies)
   set(${dependencies} ${_dependencies} PARENT_SCOPE)
 endfunction()
 
 function(_package_unset_dependencies pkg_name)
   _package_variable_unset(DEPENDENCIES ${pkg_name})
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Functions to handle reverse dependencies
 # ------------------------------------------------------------------------------
 function(_package_set_rdependencies pkg_name)
   _package_set_variable(RDEPENDENCIES ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_rdependencies pkg_name rdependencies)
   _package_get_variable(RDEPENDENCIES ${pkg_name} _rdependencies)
   set(${rdependencies} ${_rdependencies} PARENT_SCOPE)
 endfunction()
 
 function(_package_add_rdependency pkg_name rdep)
   # store the reverse dependency
   _package_add_to_variable(RDEPENDENCIES ${pkg_name} ${rdep})
 endfunction()
 
 function(_package_remove_rdependency pkg_name rdep)
   _package_remove_from_variable(RDEPENDENCIES ${pkg_name} ${rdep})
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Function to handle forcing dependencies (Package turn ON that enforce their
 # dependencies ON)
 # ------------------------------------------------------------------------------
 function(_package_set_fdependencies pkg_name)
   _package_set_variable(FDEPENDENCIES ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_fdependencies pkg_name fdependencies)
   _package_get_variable(FDEPENDENCIES ${pkg_name} _fdependencies)
   set(${fdependencies} ${_fdependencies} PARENT_SCOPE)
 endfunction()
 
 function(_package_add_fdependency pkg_name fdep)
   # store the reverse dependency
   _package_add_to_variable(FDEPENDENCIES ${pkg_name} ${fdep})
 endfunction()
 
 function(_package_remove_fdependency pkg_name fdep)
   _package_remove_from_variable(FDEPENDENCIES ${pkg_name} ${fdep})
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Documentation related functions
 # ------------------------------------------------------------------------------
 function(_package_set_documentation_files pkg_name)
   _package_set_variable(DOCUMENTATION_FILES ${pkg_name} ${ARGN})
 endfunction()
 
 function(_package_get_documentation_files pkg_name doc_files)
   _package_get_variable(DOCUMENTATION_FILES ${pkg_name} _doc_files "")
   set(${doc_files} ${_doc_files} PARENT_SCOPE)
 endfunction()
 
 function(_package_set_documentation pkg_name)
   # \n replaced by && and \\ by ££ to avoid cache problems
   set(_doc_str "")
   foreach(_str ${ARGN})
     set(_doc_str "${_doc_str}&&${_str}")
   endforeach()
 
   string(REPLACE "\\" "££" _doc_escaped "${_doc_str}")
 
   _package_set_variable(DOCUMENTATION ${pkg_name} "${_doc_str}")
 endfunction()
 
 function(_package_get_documentation pkg_name _doc)
   # \n replaced by && and \\ by ££ to avoid cache problems
   _package_get_variable(DOCUMENTATION ${pkg_name} _doc_tmp "")
 
   string(REPLACE "££" "\\" _doc_escaped "${_doc_tmp}")
   string(REPLACE "&&" "\n" _doc_newlines "${_doc_escaped}")
   set(${_doc} "${_doc_newlines}" PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Special boost thingy
 # ------------------------------------------------------------------------------
 function(_package_set_boost_component_needed pkg_name)
   _package_add_to_variable(BOOST_COMPONENTS_NEEDED ${pkg_name} ${ARGN})
   package_get_name(Boost _boost_pkg_name)
   _package_add_dependencies(${pkg_name} ${_boost_pkg_name})
 endfunction()
 
 function(_package_get_boost_component_needed pkg_name needed)
   _package_get_variable(BOOST_COMPONENTS_NEEDED ${pkg_name} _needed)
   set(${needed} ${_needed} PARENT_SCOPE)
 endfunction()
 
 function(_package_load_boost_components)
   string(TOUPPER ${PROJECT_NAME} _project)
 
   _package_get_variable_for_activated(BOOST_COMPONENTS_NEEDED _boost_needed_components)
   package_get_name(Boost _pkg_name)
 
   if(_boost_needed_components)
     message(STATUS "Looking for Boost liraries: ${_boost_needed_components}")
     foreach(_comp ${_boost_needed_components})
       find_package(Boost COMPONENTS ${_comp} QUIET)
       string(TOUPPER ${_comp} _u_comp)
       if(Boost_${_u_comp}_FOUND)
 	message(STATUS "   ${_comp}: FOUND")
 	package_set_project_variable(BOOST_${_u_comp} TRUE)
 
 	# Generate the libraries for the package
 	_package_add_libraries(${_pkg_name} ${Boost_${_u_comp}_LIBRARY})
       else()
 	message(STATUS "   ${_comp}: NOT FOUND")
       endif()
     endforeach()
   endif()
 endfunction()
 
 # ------------------------------------------------------------------------------
 # get the list of source files for a given package
 # ------------------------------------------------------------------------------
 function(_package_get_source_files pkg_name SRCS PUBLIC_HEADERS PRIVATE_HEADERS)
   string(TOUPPER ${PROJECT_NAME} _project)
 
   set(tmp_SRCS)
   set(tmp_PUBLIC_HEADERS)
   set(tmp_PRIVATE_HEADERS)
 
   foreach(_type SRCS PUBLIC_HEADERS PRIVATE_HEADERS)
     foreach(_file ${${pkg_name}_${_type}})
       string(REPLACE "${CMAKE_CURRENT_SOURCE_DIR}/" "" _rel_file "${_file}")
       list(APPEND tmp_${_type} "${_rel_file}")
     endforeach()
   endforeach()
 
   set(${SRCS}            ${tmp_SRCS}            PARENT_SCOPE)
   set(${PUBLIC_HEADERS}  ${tmp_PUBLIC_HEADERS}  PARENT_SCOPE)
   set(${PRIVATE_HEADERS} ${tmp_PRIVATE_HEADERS} PARENT_SCOPE)
 endfunction()
 
 # ==============================================================================
 # Internal functions
 # ==============================================================================
 # ------------------------------------------------------------------------------
 # Build the reverse dependencies from the dependencies
 # ------------------------------------------------------------------------------
 function(_package_build_rdependencies)
   package_get_all_packages(_pkg_list)
 
   # set empty lists
   foreach(_pkg_name ${_pkg_list})
     set(${_pkg_name}_rdeps)
   endforeach()
 
   # fill the dependencies list
   foreach(_pkg_name ${_pkg_list})
     _package_get_dependencies(${_pkg_name} _deps)
     foreach(_dep_name ${_deps})
       list(APPEND ${_dep_name}_rdeps ${_pkg_name})
     endforeach()
   endforeach()
 
   # clean and set the reverse dependencies
   foreach(_pkg_name ${_pkg_list})
     if(${_pkg_name}_rdeps)
       list(REMOVE_DUPLICATES ${_pkg_name}_rdeps)
       _package_set_rdependencies(${_pkg_name} ${${_pkg_name}_rdeps})
     endif()
   endforeach()
 endfunction()
 
 # ------------------------------------------------------------------------------
 # This function resolve the dependance order run the needed find_packages
 # ------------------------------------------------------------------------------
 function(_package_load_packages)
   package_get_all_packages(_pkg_list)
 
   # Activate the dependencies of activated package and generate an ordered list
   # of dependencies
   set(ordered_loading_list)
   foreach(_pkg_name ${_pkg_list})
     _package_load_dependencies_package(${_pkg_name} ordered_loading_list)
   endforeach()
 
   # Load the packages in the propoer order
   foreach(_pkg_name ${ordered_loading_list})
     _package_get_option_name(${_pkg_name} _option_name)
 
     if(${_option_name})
       _package_load_package(${_pkg_name})
     else()
       # deactivate the packages than can already be deactivated
       _package_deactivate(${_pkg_name})
     endif()
   endforeach()
 
   # generates the activated and unactivated lists of packages
   set(_packages_activated)
   set(_packages_deactivated)
   foreach(_pkg_name ${_pkg_list})
     _package_is_activated(${_pkg_name} _active)
     if(_active)
       list(APPEND _packages_activated ${_pkg_name})
     else()
       list(APPEND _packages_deactivated ${_pkg_name})
     endif()
   endforeach()
 
   # generate the list usable by the calling code
   package_set_project_variable(ACTIVATED_PACKAGE_LIST "${_packages_activated}")
   package_set_project_variable(DEACTIVATED_PACKAGE_LIST "${_packages_deactivated}")
 endfunction()
 
 
 
 # ------------------------------------------------------------------------------
 # This load an external package and recursively all its dependencies
 # ------------------------------------------------------------------------------
 function(_package_load_dependencies_package pkg_name loading_list)
   # Get packages informations
   _package_get_option_name(${pkg_name} _pkg_option_name)
   _package_get_dependencies(${pkg_name} _dependencies)
 
   # handle the dependencies
   foreach(_dep_name ${_dependencies})
     _package_get_description(${_dep_name} _dep_desc)
     _package_get_option_name(${_dep_name} _dep_option_name)
 
     _package_get_fdependencies(${_dep_name} _fdeps)
     if(${_pkg_option_name})
       if("${_fdeps}" STREQUAL "")
 	set(${_dep_name}_OLD ${${_dep_option_name}} CACHE INTERNAL "" FORCE)
       endif()
 
       # set the option to on
       set(${_dep_option_name} ON CACHE BOOL "${_dep_desc}" FORCE)
 
       # store the reverse dependency
       _package_add_fdependency(${_dep_name} ${pkg_name})
     else()
       # check if this is the last reverse dependency
       list(LENGTH _fdeps len)
       list(FIND _fdeps ${pkg_name} pos)
       if((len EQUAL 1) AND (NOT pos EQUAL -1))
 	set(${_dep_option_name} ${${_dep_name}_OLD} CACHE BOOL "${_dep_desc}" FORCE)
 	unset(${_dep_name}_OLD CACHE)
       endif()
 
       # remove the pkg_name form the reverse dependency
       _package_remove_fdependency(${_dep_name} ${pkg_name})
     endif()
 
     # recusively load the dependencies
     _package_load_dependencies_package(${_dep_name} ${loading_list})
   endforeach()
 
   # get the compile flags
   _package_get_compile_flags(${pkg_name} _pkg_comile_flags)
 
   # if package option is on add it in the list
   if(${_pkg_option_name})
     list(FIND ${loading_list} ${pkg_name} _pos)
     if(_pos EQUAL -1)
       set(_tmp_loading_list ${${loading_list}})
       list(APPEND _tmp_loading_list ${pkg_name})
       set(${loading_list} "${_tmp_loading_list}" PARENT_SCOPE)
     endif()
 
     #add the comilation flags if needed
     if(_pkg_comile_flags)
       add_flags(cxx ${_pkg_comile_flags})
     endif()
   else()
     #remove the comilation flags if needed
     if(_pkg_comile_flags)
       remove_flags(cxx ${_pkg_comile_flags})
     endif()
   endif()
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Load the package if it is an external one
 # ------------------------------------------------------------------------------
 function(_package_load_package pkg_name)
   # load the package if it is an external
   _package_get_nature(${pkg_name} _nature)
   if(${_nature} MATCHES "external")
     _package_use_system(${pkg_name} _use_system)
 
     set(_activated TRUE)
     if(_use_system)
       _package_load_external_package(${pkg_name} _activated)
     else()
       _package_load_third_party_script(${pkg_name})
 
       string(TOUPPER ${${pkg_name}} _u_package)
       if(${_u_package}_LIBRARIES)
 	_package_set_libraries(${pkg_name} ${${_u_package}_LIBRARIES})
       endif()
       if(${_u_package}_INCLUDE_DIR)
 	_package_set_include_dir(${pkg_name} ${${_u_package}_INCLUDE_DIR})
       endif()
     endif()
 
     if(_activated)
       _package_activate(${pkg_name})
     elseif()
       _package_deactivate(${pkg_name})
     endif()
   else(${_nature})
     _package_activate(${pkg_name})
   endif()
 endfunction()
 
 
 # ------------------------------------------------------------------------------
 # Load external packages
 # ------------------------------------------------------------------------------
 function(_package_load_external_package pkg_name activate)
   _package_get_find_package_extra_options(${pkg_name} _options)
   if(_options)
     cmake_parse_arguments(_opt_pkg "" "LANGUAGE" "PREFIX;FOUND;ARGS" ${_options})
     if(_opt_pkg_UNPARSED_ARGUMENTS)
       message("You passed too many options for the find_package related to ${${pkg_name}} \"${_opt_pkg_UNPARSED_ARGUMENTS}\"")
     endif()
   endif()
 
   if(_opt_pkg_LANGUAGE)
     foreach(_language ${_opt_pkg_LANGUAGE})
       enable_language(${_language})
     endforeach()
   endif()
 
   _package_get_real_name(${pkg_name} _real_name)
 
   # find the package
   find_package(${_real_name} REQUIRED ${_opt_pkg_ARGS})
 
   # check if the package is found
   if(_opt_pkg_PREFIX)
     set(_package_prefix ${_opt_pkg_PREFIX})
   else()
     string(TOUPPER ${${pkg_name}} _u_package)
     set(_package_prefix ${_u_package})
   endif()
 
   set(_act FALSE)
   set(_prefix_to_consider)
-  if(_opt_pkg_FOUND)
-    set(_act TRUE)
-    set(_prefix_to_consider ${_package_prefix})
+  if(DEFINED _opt_pkg_FOUND)
+    if(${_opt_pkg_FOUND})
+      set(_act TRUE)
+      set(_prefix_to_consider ${_package_prefix})
+    endif()
   else()
     foreach(_prefix ${_package_prefix})
       if(${_prefix}_FOUND)
 	set(_act TRUE)
 	list(APPEND _prefix_to_consider ${_prefix})
       endif()
     endforeach()
   endif()
 
   if(_act)
     foreach(_prefix ${_prefix_to_consider})
       # Generate the include dir for the package
       if(DEFINED ${_prefix}_INCLUDE_DIRS)
 	_package_set_include_dir(${_pkg_name} ${${_prefix}_INCLUDE_DIRS})
       elseif(DEFINED ${_prefix}_INCLUDE_DIR)
 	_package_set_include_dir(${_pkg_name} ${${_prefix}_INCLUDE_DIR})
       elseif(DEFINED ${_prefix}_INCLUDE_PATH)
 	_package_set_include_dir(${_pkg_name} ${${_prefix}_INCLUDE_PATH})
       endif()
 
       # Generate the libraries for the package
       if(DEFINED ${_prefix}_LIBRARIES)
 	_package_set_libraries(${_pkg_name} ${${_prefix}_LIBRARIES})
       elseif(DEFINED ${_prefix}_LIBRARY)
 	_package_set_libraries(${_pkg_name} ${${_prefix}_LIBRARY})
       endif()
     endforeach()
   endif()
   set(${activate} ${_act} PARENT_SCOPE)
 endfunction()
 
 # ------------------------------------------------------------------------------
 # Sanity check functions
 # ------------------------------------------------------------------------------
 function(_package_check_files_exists)
   set(_message FALSE)
 
   package_get_all_packages(_pkg_list)
   foreach(_pkg_name ${_pkg_list})
     set(_pkg_files
       ${${_pkg_name}_SRCS}
       ${${_pkg_name}_PUBLIC_HEADERS}
       ${${_pkg_name}_PRIVATE_HEADERS}
       )
 
     _package_get_real_name(${_pkg_name} _real_name)
 
     foreach(_file ${_pkg_files})
       if(NOT EXISTS "${_file}")
 	if(NOT _message)
 	  set(_message TRUE)
 	  message("This file(s) is(are) present in a package but are not present on disk.")
 	endif()
 
 	message(" PACKAGE ${_real_name} FILE ${_file}")
       endif()
     endforeach()
   endforeach()
 
   if(_message)
     message(SEND_ERROR "Please check the content of your packages to correct this warnings")
   endif()
 endfunction()
 
 # ------------------------------------------------------------------------------
 function(_package_check_files_registered)
   set(_pkg_files)
   package_get_all_packages(_pkg_list)
   # generates a file list of registered files
   foreach(_pkg_name ${_pkg_list})
     list(APPEND _pkg_files
       ${${_pkg_name}_SRCS}
       ${${_pkg_name}_PUBLIC_HEADERS}
       ${${_pkg_name}_PRIVATE_HEADERS}
       )
   endforeach()
 
   # generates the list of files in the source folders
   set(_all_src_files)
   foreach(_src_folder ${ARGN})
     foreach(_ext "cc" "hh" "c" "h" "hpp")
       file(GLOB_RECURSE _src_files "${_src_folder}/*.${_ext}")
       list(APPEND _all_src_files ${_src_files})
     endforeach()
   endforeach()
 
   if(_all_src_files)
     list(REMOVE_DUPLICATES _all_src_files)
   endif()
 
   set(_not_registerd_files)
   # check only sources files ine the source folders
   foreach(_src_folder ${ARGN})
     foreach(_file ${_all_src_files})
       if("${_file}" MATCHES  "${_src_folder}")
 	list(FIND _pkg_files "${_file}" _index)
 	if (_index EQUAL -1)
 	  list(APPEND _not_registerd_files ${_file})
 	endif()
       endif()
     endforeach()
   endforeach()
 
   if(AUTO_MOVE_UNKNOWN_FILES)
     file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/unknown_files)
   endif()
 
   # warn the user and move the files if needed
   if(_not_registerd_files)
     if(EXISTS ${PROJECT_BINARY_DIR}/missing_files_in_packages)
       file(REMOVE ${PROJECT_BINARY_DIR}/missing_files_in_packages)
     endif()
 
     message("This files are present in the source folders but are not registered in any package")
     foreach(_file ${_not_registerd_files})
       message(" ${_file}")
       if(AUTO_MOVE_UNKNOWN_FILES)
 	get_filename_component(_file_name ${_file} NAME)
 	file(RENAME ${_file} ${PROJECT_BINARY_DIR}/unknown_files/${_file_name})
       endif()
 
       file(APPEND ${PROJECT_BINARY_DIR}/missing_files_in_packages "${_file}
 ")
     endforeach()
 
     if(AUTO_MOVE_UNKNOWN_FILES)
       message(SEND_ERROR "The files where moved in the followinf folder ${PROJECT_BINARY_DIR}/unknown_files\n
 Please register them in the good package or clean the sources")
     else()
       message(SEND_ERROR "Please register them in the good package or clean the sources")
     endif()
 
   endif()
 endfunction()
 
 # ------------------------------------------------------------------------------
\ No newline at end of file
diff --git a/doc/manual/manual-cohesive_elements_insertion.tex b/doc/manual/manual-cohesive_elements_insertion.tex
index 31a50c17c..a13692960 100644
--- a/doc/manual/manual-cohesive_elements_insertion.tex
+++ b/doc/manual/manual-cohesive_elements_insertion.tex
@@ -1,85 +1,81 @@
 For cohesive material, \akantu has a pre-defined material selector to assign 
 the first cohesive material by default to the cohesive elements which is called 
 \code{DefaultMaterialCohesiveSelector} and it inherits its properties from 
 \code{DefaultMaterialSelector}. Multiple cohesive materials can be assigned 
 using mesh data information (for more details, see \ref{intrinsic_insertion}).
 
 \subsection{Insertion of Cohesive Elements}
+<<<<<<< Updated upstream
 Cohesive elements are currently compatible only with static simulation
 and dynamic simulation with an explicit time integration scheme (see
-section~\ref{ssect:smm:expl-time-integr}). They do not have to be
-inserted when the mesh is generated but during the simulation. At any
-time during the simulation, it is possible to access the following
-energies with the relative function:
+section~\ref{ssect:smm:expl-time-integr}). They do not have to be inserted when the mesh is generated (intrinsic) but can be added during the simulation (extrinsic). At any time during the simulation, it is possible to access the following energies with the relative function:
 \begin{cpp}
   Real Ed = model.getEnergy("dissipated");
   Real Er = model.getEnergy("reversible");
   Real Ec = model.getEnergy("contact");
 \end{cpp}
 
+A new model have to be call in a very similar way that the solid mechanics model:
+\begin{cpp}
+  SolidMechanicsModelCohesive model(mesh);
+  model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
+\end{cpp} 
+
+
 \subsubsection{Extrinsic approach}
 The dynamic insertion of extrinsic cohesive elements should be initialized 
 in the following way:
 \begin{cpp}
-  SolidMechanicsModelCohesive model(mesh);
-  model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
   model.updateAutomaticInsertion();
 \end{cpp} 
-During the simulation, stress has to be checked along each facet in order to 
-eventually insert cohesive elements.
+During the simulation, stress has to be checked along each facet in order to insert cohesive elements where thestress criterion is reached.
 This check is performed by calling the method \code{checkCohesiveStress}, as 
 example before each step resolution:
 \begin{cpp}
   model.checkCohesiveStress();
   model.solveStep();
 \end{cpp}
 The area where stresses are checked and cohesive elements inserted can be limited 
 using the method \code{limitInsertion} during initialization. As example, to 
 limit insertion in the range $[-1.5, 1.5]$ in the $x$ direction: 
 \begin{cpp}
-  SolidMechanicsModelCohesive model(mesh);
-  model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
   model.limitInsertion(_x, -1.5, 1.5);
   model.updateAutomaticInsertion();
 \end{cpp} 
 Additional restrictions with respect to $y$ and $z$ directions can be added as well.
 
 \subsubsection{Intrinsic approach \label{intrinsic_insertion}}
 Intrinsic cohesive elements are inserted in the mesh with the method 
 \code{insertIntrinsicElements}. Similarly, the range of insertion can me limited 
 with \code{limitInsertion}. As example with a static simulation,
 \begin{cpp}
-  SolidMechanicsModelCohesive model(mesh);
-  model.initFull(SolidMechanicsModelCohesiveOptions(_static));
   model.limitInsertion(_x, -1.5, 1.5);
   model.insertIntrinsicElements();
 \end{cpp} 
 Mesh data information becomes vital to the insertion of cohesive elements along 
 surface with more sophisticated geometry or when multiple cohesive materials are 
-wanted. At this purpose, cohesive elements can be easily inserted along a 
+wanted. To do so, cohesive elements can be inserted along a 
 specific group of surface elements identified in a GMSH geometry file. This can 
-be achieved, in the input file, by specifying the name of these physical groups 
-in the \textit{mesh parameters} section as well as their corresponding cohesive 
-materials. As example, with two physical surfaces named 
+be achieved with the material selector (see section~\ref{sect:smm:materialselector}), in the input file specify the name of these physical groups in the corresponding cohesive materials, and call these material in the \textit{mesh parameters} section. As example, with two physical surfaces named 
 \textit{weak\_interface} and \textit{strong\_interface} defined in the GMSH 
 geometry file:
 \begin{cpp}
 ...
   material %\emph{cohesive constitutive\_law}% [
      name = weak_interface
      sigma_c = $value$
      ...
   ]
   material %\emph{cohesive constitutive\_law}% [
      name = strong_interface
      sigma_c = $value$
      ...
   ]
   mesh parameters [
      	cohesive_surfaces = weak_interface,strong_interface
   ]
 \end{cpp}
 
 In this case, there is no need to call \code{insertIntrinsicElements} anymore 
 since the insertion of cohesive elements along physical surfaces is performed 
 automatically during \code{initFull} call.    
diff --git a/packages/blackdynamite.cmake b/packages/blackdynamite.cmake
index 5b973c202..b4ef8ee38 100644
--- a/packages/blackdynamite.cmake
+++ b/packages/blackdynamite.cmake
@@ -1,41 +1,41 @@
 #===============================================================================
 # @file   blackdynamite.cmake
 #
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @date   Tue Nov 29 15:16:35 2011
 #
 # @brief  package description for BlackDynamite support
 #
 # @section LICENSE
 #
 # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # Akantu is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Akantu. If not, see <http://www.gnu.org/licenses/>.
 #
 #===============================================================================
 package_declare(BlackDynamite EXTERNAL
   DESCRIPTION "Use BlackDynamite library"
-  SYSTEM OFF)
+  SYSTEM OFF third-party/cmake/blackdynamite.cmake
+  EXTRA_PACKAGE_OPTIONS FOUND BlackDynamite_FOUND)
 
-
-set(_version 1)
+set(_version master)
 
 package_add_third_party_script_variable(BlackDynamite
   BLACKDYNAMITE_VERSION "${_version}")
 package_add_third_party_script_variable(BlackDynamite
-  BLACKDYNAMITE_GIT "svn+ssh://lsmssrv1.epfl.ch/space/repositories/SimulPack/BlackDynamite")
+  BLACKDYNAMITE_GIT "git@lsmssrv1.epfl.ch:blackdynamite.git")
 package_add_third_party_script_variable(BlackDynamite
-  BLACKDYNAMITE_ARCHIVE "blackdynamite-${_version}.tar.gz")
\ No newline at end of file
+  BLACKDYNAMITE_ARCHIVE "blackdynamite-${_version}.tar.gz")
diff --git a/src/geometry/mesh_abstract_intersector.hh b/src/geometry/mesh_abstract_intersector.hh
index 17060ee9c..ad7dd115d 100644
--- a/src/geometry/mesh_abstract_intersector.hh
+++ b/src/geometry/mesh_abstract_intersector.hh
@@ -1,78 +1,106 @@
 /**
  * @file mesh_abstract_intersector.hh
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date creation: Mon Jul 13 2015
  * @date last modification: Mon Jul 13 2015
  *
  * @brief Abstract class for intersection computations
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_ABSTRACT_INTERSECTOR_HH__
 #define __AKANTU_MESH_ABSTRACT_INTERSECTOR_HH__
 
 #include "aka_common.hh"
 #include "mesh_geom_abstract.hh"
 
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_AKANTU__
 
 /**
  * @brief Class used to perform intersections on a mesh and construct output data
  */
 template<class Query>
 class MeshAbstractIntersector : public MeshGeomAbstract {
 
 public:
   /// Construct from mesh
   explicit MeshAbstractIntersector(Mesh & mesh);
 
   /// Destructor
   virtual ~MeshAbstractIntersector();
 
 public:
+  /* ------------------------------------------------------------------------ */
+  /* Accessors                                                                */
+  /* ------------------------------------------------------------------------ */
+public:
+  /// get the new_node_per_elem array
+  AKANTU_GET_MACRO(NewNodePerElem, new_node_per_elem, const Array<UInt> *)
+  AKANTU_GET_MACRO_NOT_CONST(NewNodePerElemNotConst, new_node_per_elem, Array<UInt> *)
+  /// get the new_nodes array
+  AKANTU_GET_MACRO(NewNodes, new_nodes, const Array<Real> *)
+  /// get the nb_seg_by_el UInt
+  AKANTU_GET_MACRO(NbSegByEl, nb_seg_by_el, const UInt)
+
   /**
    * @brief Compute the intersection with a query object
    *
    * This function needs to be implemented for every subclass. It computes the intersections
    * with the tree of primitives and creates the data for the user.
    * 
    * @param query the CGAL primitive of the query object
    */
   virtual void computeIntersectionQuery(const Query & query) = 0;
 
-  /// Compute list of queries
+  /// Compute intersection points between the mesh and a query (closed surface in 3D or a closed curve in 2D)
+  virtual void computeMeshQueryIntersectionPoint(const Query & query) = 0;
+
+  /// Compute intersection between the mesh and a list of queries
   virtual void computeIntersectionQueryList(const std::list<Query> & query_list);
 
+  /// Compute intersection points between the mesh and a list of queries
+  virtual void  computeMeshQueryListIntersectionPoint(const std::list<Query> & query_list);
+
   /// Compute whatever result is needed from the user
   virtual void buildResultFromQueryList(const std::list<Query> & query_list) = 0;
+
+protected:
+  /// new node per element (column 0: number of new nodes, then odd is the intersection node number and even the ID of the intersected segment)
+  Array<UInt> * new_node_per_elem;
+
+  /// intersection output: new intersection points (computeMeshQueryListIntersectionPoint)
+  Array<Real> * new_nodes;
+
+  /// number of segment in a considered element of the templated type of element specialized intersector
+  const UInt nb_seg_by_el;
 };
 
 __END_AKANTU__
 
 #include "mesh_abstract_intersector_tmpl.hh"
 
 #endif // __AKANTU_MESH_ABSTRACT_INTERSECTOR_HH__
diff --git a/src/geometry/mesh_abstract_intersector_tmpl.hh b/src/geometry/mesh_abstract_intersector_tmpl.hh
index 3d809e776..951ff2874 100644
--- a/src/geometry/mesh_abstract_intersector_tmpl.hh
+++ b/src/geometry/mesh_abstract_intersector_tmpl.hh
@@ -1,68 +1,87 @@
 /**
  * @file mesh_abstract_intersector_tmpl.hh
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date creation: Mon Jul 13 2015
  * @date last modification: Mon Jul 13 2015
  *
  * @brief General class for intersection computations
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_ABSTRACT_INTERSECTOR_TMPL_HH__
 #define __AKANTU_MESH_ABSTRACT_INTERSECTOR_TMPL_HH__
 
 #include "aka_common.hh"
 #include "mesh_abstract_intersector.hh"
 
 __BEGIN_AKANTU__
 
 template<class Query>
 MeshAbstractIntersector<Query>::MeshAbstractIntersector(Mesh & mesh):
-  MeshGeomAbstract(mesh)
+  MeshGeomAbstract(mesh),
+  new_node_per_elem(NULL),
+  new_nodes(NULL),
+  nb_seg_by_el(0)
 {}
 
 template<class Query>
 MeshAbstractIntersector<Query>::~MeshAbstractIntersector()
 {}
 
 template<class Query>
 void MeshAbstractIntersector<Query>::computeIntersectionQueryList(
   const std::list<Query> & query_list) {
   AKANTU_DEBUG_IN();
   
   typename std::list<Query>::const_iterator
     query_it = query_list.begin(),
     query_end = query_list.end();
 
   for (; query_it != query_end ; ++query_it) {
     computeIntersectionQuery(*query_it);
   }
   
   AKANTU_DEBUG_OUT();
 }
 
+template<class Query>
+void MeshAbstractIntersector<Query>::computeMeshQueryListIntersectionPoint(
+  const std::list<Query> & query_list) {
+  AKANTU_DEBUG_IN();
+  
+  typename std::list<Query>::const_iterator
+    query_it = query_list.begin(),
+    query_end = query_list.end();
+
+  for (; query_it != query_end ; ++query_it) {
+    computeMeshQueryIntersectionPoint(*query_it);
+  }
+  
+  AKANTU_DEBUG_OUT();
+}
+
 __END_AKANTU__
 
 #endif // __AKANTU_MESH_ABSTRACT_INTERSECTOR_TMPL_HH__
diff --git a/src/geometry/mesh_geom_intersector_tmpl.hh b/src/geometry/mesh_geom_intersector_tmpl.hh
index 875aa366c..4e33228a7 100644
--- a/src/geometry/mesh_geom_intersector_tmpl.hh
+++ b/src/geometry/mesh_geom_intersector_tmpl.hh
@@ -1,59 +1,62 @@
 /**
  * @file mesh_geom_intersector_tmpl.hh
  *
  * @author Lucas Frerot <lucas.frerot@epfl.ch>
  *
  * @date creation: Wed Apr 29 2015
  * @date last modification: Wed Apr 29 2015
  *
  * @brief General class for intersection computations
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_GEOM_INTERSECTOR_TMPL_HH__
 #define __AKANTU_MESH_GEOM_INTERSECTOR_TMPL_HH__
 
 #include "aka_common.hh"
 #include "mesh_geom_intersector.hh"
 
+/* -------------------------------------------------------------------------- */
+
 __BEGIN_AKANTU__
 
 template<UInt dim, ElementType type, class Primitive, class Query, class Kernel>
 MeshGeomIntersector<dim, type, Primitive, Query, Kernel>::MeshGeomIntersector(Mesh & mesh) :
   MeshAbstractIntersector<Query>(mesh),
   factory(mesh)
 {}
 
 template<UInt dim, ElementType type, class Primitive, class Query, class Kernel>
 MeshGeomIntersector<dim, type, Primitive, Query, Kernel>::~MeshGeomIntersector()
 {}
 
 template<UInt dim, ElementType type, class Primitive, class Query, class Kernel>
 void MeshGeomIntersector<dim, type, Primitive, Query, Kernel>::constructData() {
+  this->new_nodes->resize(0);
   factory.constructData();
 }
 
 __END_AKANTU__
 
 #endif // __AKANTU_MESH_GEOM_INTERSECTOR_TMPL_HH__
 
diff --git a/src/geometry/mesh_segment_intersector.hh b/src/geometry/mesh_segment_intersector.hh
index 2ecfe0f54..f857f08ba 100644
--- a/src/geometry/mesh_segment_intersector.hh
+++ b/src/geometry/mesh_segment_intersector.hh
@@ -1,98 +1,101 @@
 /**
  * @file mesh_segment_intersector.hh
  *
  * @author Lucas Frerot
  *
  * @date creation: Wed Apr 29 2015
  * @date last modification: Wed Apr 29 2015
  *
  * @brief Computation of mesh intersection with segments
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_SEGMENT_INTERSECTOR_HH__
 #define __AKANTU_MESH_SEGMENT_INTERSECTOR_HH__
 
 #include "aka_common.hh"
 #include "mesh_geom_intersector.hh"
 
 #include "mesh_geom_common.hh"
 
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_AKANTU__
 
 /// Here, we know what kernel we have to use
 typedef Cartesian K;
 
 template<UInt dim, ElementType type>
 class MeshSegmentIntersector : public MeshGeomIntersector<dim, type, Triangle<K>, K::Segment_3, K> {
   /// Parent class type
   typedef MeshGeomIntersector<dim, type, Triangle<K>, K::Segment_3, K> parent_type;
 
   /// Result of intersection function type
   typedef typename IntersectionTypeHelper<TreeTypeHelper< Triangle<K>, K>, K::Segment_3>::intersection_type result_type;
 
   /// Pair of segments and element id
   typedef std::pair<K::Segment_3, UInt> pair_type;
 
 public:
   /// Construct from mesh
   explicit MeshSegmentIntersector(Mesh & mesh, Mesh & result_mesh);
 
   /// Destructor
   virtual ~MeshSegmentIntersector();
 
 public:
   /**
    * @brief Computes the intersection of the mesh with a segment
    *
    * @param query the segment to compute the intersections with the mesh
    */
   virtual void computeIntersectionQuery(const K::Segment_3 & query);
 
+  /// Compute intersection points between the mesh and a query
+  virtual void computeMeshQueryIntersectionPoint(const K::Segment_3 & query);
+
   /// Compute the embedded mesh
   virtual void buildResultFromQueryList(const std::list<K::Segment_3> & query_list);
 
   void setPhysicalName(const std::string & other)
   { current_physical_name = other; }
 
 protected:
   /// Compute segments from intersection list
   void computeSegments(const std::list<result_type> & intersections,
                        std::set<pair_type, segmentPairsLess> & segments,
                        const K::Segment_3 & query);
 
 protected:
   /// Result mesh
   Mesh & result_mesh;
 
   /// Physical name of the current batch of queries
   std::string current_physical_name;
 };
  
 __END_AKANTU__
 
 #include "mesh_segment_intersector_tmpl.hh"
 
 #endif // __AKANTU_MESH_SEGMENT_INTERSECTOR_HH__
diff --git a/src/geometry/mesh_segment_intersector_tmpl.hh b/src/geometry/mesh_segment_intersector_tmpl.hh
index 7989a25f7..64cc4c124 100644
--- a/src/geometry/mesh_segment_intersector_tmpl.hh
+++ b/src/geometry/mesh_segment_intersector_tmpl.hh
@@ -1,261 +1,267 @@
 /**
  * @file mesh_segment_intersector_tmpl.hh
  *
  * @author Lucas Frerot
  *
  * @date creation: Wed Apr 29 2015
  * @date last modification: Wed Apr 29 2015
  *
  * @brief Computation of mesh intersection with segments
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__
 #define __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__
 
 #include "aka_common.hh"
 #include "mesh_geom_common.hh"
 #include "tree_type_helper.hh"
 
 __BEGIN_AKANTU__
 
 template<UInt dim, ElementType type>
 MeshSegmentIntersector<dim, type>::MeshSegmentIntersector(Mesh & mesh, Mesh & result_mesh):
   parent_type(mesh),
   result_mesh(result_mesh),
   current_physical_name()
 {
+  this->new_nodes = new Array<Real>(0,dim);
   this->constructData();
 }
 
 template<UInt dim, ElementType type>
 MeshSegmentIntersector<dim, type>::~MeshSegmentIntersector()
 {}
 
 template<UInt dim, ElementType type>
 void MeshSegmentIntersector<dim, type>::computeIntersectionQuery(const K::Segment_3 & query) {
   AKANTU_DEBUG_IN();
   
   result_mesh.addConnectivityType(_segment_2, _not_ghost);
   result_mesh.addConnectivityType(_segment_2, _ghost);
 
   std::list<result_type> result_list;
   std::set<std::pair<K::Segment_3, UInt>, segmentPairsLess> segment_set;
 
   this->factory.getTree().all_intersections(query, std::back_inserter(result_list));
   this->computeSegments(result_list, segment_set, query);
 
   // Arrays for storing nodes and connectivity
   Array<Real> & nodes = result_mesh.getNodes();
   Array<UInt> & connectivity = result_mesh.getConnectivity(_segment_2);
 
   // Arrays for storing associated element and physical name
   bool valid_elemental_data = true;
   Array<Element> * associated_element = NULL;
   Array<std::string> * associated_physical_name = NULL;
 
   try {
     associated_element = &result_mesh.getData<Element>("associated_element", _segment_2);
     associated_physical_name = &result_mesh.getData<std::string>("physical_names", _segment_2);
   } catch (debug::Exception & e) {
     valid_elemental_data = false;
   }
 
   std::set<pair_type, segmentPairsLess>::iterator
     it = segment_set.begin(),
     end = segment_set.end();
 
   // Loop over the segment pairs
   for (; it != end ; ++it) {
     if (!it->first.is_degenerate()) {
       Vector<UInt> segment_connectivity(2);
       segment_connectivity(0) = result_mesh.getNbNodes();
       segment_connectivity(1) = result_mesh.getNbNodes() + 1;
       connectivity.push_back(segment_connectivity);
 
       // Copy nodes
       Vector<Real> source(dim), target(dim);
       for (UInt j = 0 ; j < dim ; j++) {
         source(j) = it->first.source()[j];
         target(j) = it->first.target()[j];
       }
 
       nodes.push_back(source);
       nodes.push_back(target);
 
       // Copy associated element info
       if (valid_elemental_data) {
         associated_element->push_back(Element(type, it->second));
         associated_physical_name->push_back(current_physical_name);
       }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
+template<UInt dim, ElementType type>
+void MeshSegmentIntersector<dim, type>:: computeMeshQueryIntersectionPoint(const K::Segment_3 & query) {
+  AKANTU_DEBUG_ERROR("The method: computeMeshQueryIntersectionPoint has not been implemented in class MeshSegmentIntersector!");
+}
+
 template<UInt dim, ElementType type>
 void MeshSegmentIntersector<dim, type>::buildResultFromQueryList(const std::list<K::Segment_3> & query_list) {
   AKANTU_DEBUG_IN();
 
   this->computeIntersectionQueryList(query_list);
 
   AKANTU_DEBUG_OUT();
 }
 
 template<UInt dim, ElementType type>
 void MeshSegmentIntersector<dim, type>::computeSegments(const std::list<result_type> & intersections,
                                                         std::set<pair_type, segmentPairsLess> & segments,
                                                         const K::Segment_3 & query) {
   AKANTU_DEBUG_IN();
   
   /*
    * Number of intersections = 0 means
    *
    * - query is completely outside mesh
    * - query is completely inside primitive
    *
    * We try to determine the case and still construct the segment list
    */
   if (intersections.size() == 0) {
     // We look at all the primitives intersected by two rays
     // If there is one primitive in common, then query is inside
     // that primitive
     K::Ray_3 ray1(query.source(), query.target());
     K::Ray_3 ray2(query.target(), query.source());
 
     std::set<UInt> ray1_results, ray2_results;
 
     this->factory.getTree().all_intersected_primitives(ray1, std::inserter(ray1_results, ray1_results.begin()));
     this->factory.getTree().all_intersected_primitives(ray2, std::inserter(ray2_results, ray2_results.begin()));
 
     bool inside_primitive = false;
     UInt primitive_id = 0;
 
     std::set<UInt>::iterator
       ray2_it  = ray2_results.begin(),
       ray2_end = ray2_results.end();
 
     // Test if first list contains an element of second list
     for (; ray2_it != ray2_end && !inside_primitive ; ++ray2_it) {
       if (ray1_results.find(*ray2_it) != ray1_results.end()) {
         inside_primitive = true;
         primitive_id = *ray2_it;
       }
     }
 
     if (inside_primitive) {
       segments.insert(std::make_pair(query, primitive_id));
     }
   }
 
   else {
     typename std::list<result_type>::const_iterator
       it = intersections.begin(),
       end = intersections.end();
 
     for(; it != end ; ++it) {
       UInt el = (*it)->second;
 
       // Result of intersection is a segment
       if (const K::Segment_3 * segment = boost::get<K::Segment_3>(&((*it)->first))) {
         // Check if the segment was alread created
         segments.insert(std::make_pair(*segment, el));
       }
 
       // Result of intersection is a point
       else if (const K::Point_3 * point = boost::get<K::Point_3>(&((*it)->first))) {
         // We only want to treat points differently if we're in 3D with Tetra4 elements
         // This should be optimized by compilator
         if (dim == 3 && type == _tetrahedron_4) {
           UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
           TreeTypeHelper<Triangle<K>, K>::container_type facets;
 
           const Array<Real> & nodes = this->mesh.getNodes();
           Array<UInt>::const_vector_iterator
             connectivity_vec = this->mesh.getConnectivity(type).begin(nb_nodes_per_element);
 
           const Vector<UInt> & el_connectivity = connectivity_vec[el];
 
           Matrix<Real> node_coordinates(dim, nb_nodes_per_element);
           for (UInt i = 0 ; i < nb_nodes_per_element ; i++)
             for (UInt j = 0 ; j < dim ; j++)
               node_coordinates(j, i) = nodes(el_connectivity(i), j);
 
           this->factory.addPrimitive(node_coordinates, el, facets);
 
           // Local tree
           TreeTypeHelper<Triangle<K>, K>::tree * local_tree =
             new TreeTypeHelper<Triangle<K>, K>::tree(facets.begin(), facets.end());
 
           // Compute local intersections (with current element)
           std::list<result_type> local_intersections;
           local_tree->all_intersections(query, std::back_inserter(local_intersections));
 
           bool out_point_found = false;
           typename std::list<result_type>::const_iterator
             local_it = local_intersections.begin(),
             local_end = local_intersections.end();
 
           for (; local_it != local_end ; ++local_it) {
             if (const K::Point_3 * local_point = boost::get<K::Point_3>(&((*local_it)->first))) {
               if (!comparePoints(*point, *local_point)) {
                 K::Segment_3 seg(*point, *local_point);
                 segments.insert(std::make_pair(seg, el));
                 out_point_found = true;
               }
             }
           }
 
           if (!out_point_found) {
             TreeTypeHelper<Triangle<K>, K>::point_type
               a(node_coordinates(0, 0), node_coordinates(1, 0), node_coordinates(2, 0)),
               b(node_coordinates(0, 1), node_coordinates(1, 1), node_coordinates(2, 1)),
               c(node_coordinates(0, 2), node_coordinates(1, 2), node_coordinates(2, 2)),
               d(node_coordinates(0, 3), node_coordinates(1, 3), node_coordinates(2, 3));
             K::Tetrahedron_3 tetra(a, b, c, d);
             const K::Point_3 * inside_point = NULL;
             if (tetra.has_on_bounded_side(query.source()) && !tetra.has_on_boundary(query.source()))
               inside_point = &query.source();
             else if (tetra.has_on_bounded_side(query.target()) && !tetra.has_on_boundary(query.target()))
               inside_point = &query.target();
 
             if (inside_point) {
               K::Segment_3 seg(*inside_point, *point);
               segments.insert(std::make_pair(seg, el));
             }
           }
 
           delete local_tree;
         }
       }
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 __END_AKANTU__
 
 #endif // __AKANTU_MESH_SEGMENT_INTERSECTOR_TMPL_HH__
 
diff --git a/src/geometry/mesh_sphere_intersector.hh b/src/geometry/mesh_sphere_intersector.hh
index 55bd08ab7..40a0651b8 100644
--- a/src/geometry/mesh_sphere_intersector.hh
+++ b/src/geometry/mesh_sphere_intersector.hh
@@ -1,119 +1,118 @@
 /**
- * @file mesh_segment_intersector.hh
+ * @file mesh_sphere_intersector.hh
  *
  * @author Clement Roux-Langlois <clement.roux@epfl.ch>
  *
  * @date creation: Wed Jun 10 2015
  *
  * @brief Computation of mesh intersection with sphere(s)
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_HH__
 #define __AKANTU_MESH_SPHERE_INTERSECTOR_HH__
 
 #include "aka_common.hh"
 #include "mesh_geom_intersector.hh"
 
 #include "mesh_geom_common.hh"
 
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_AKANTU__
 
 /* -------------------------------------------------------------------------- */
 /* class for new igfem elements mesh events                                   */
 /* -------------------------------------------------------------------------- */
 #if defined(AKANTU_IGFEM)
 class NewIGFEMElementsEvent : public NewElementsEvent {
 public:
   AKANTU_GET_MACRO_NOT_CONST(OldElementsList, old_elements, Array<Element> &);
   AKANTU_GET_MACRO(OldElementsList, old_elements, const Array<Element> &);
 protected:
   Array<Element> old_elements;
 };
 #endif
 
 
 /// Here, we know what kernel we have to use
 typedef Spherical SK;
 
 template<UInt dim, ElementType type>
 class MeshSphereIntersector : public MeshGeomIntersector<dim, type, Line_arc<SK>, SK::Sphere_3, SK> {
-  /// Parent class type
-  typedef MeshGeomIntersector<dim, type, Line_arc<SK>, SK::Sphere_3, SK> parent_type;
-
   /// Result of intersection function type
   typedef typename IntersectionTypeHelper<TreeTypeHelper< Triangle<K>, K>, K::Segment_3>::intersection_type result_type;
 
   /// Pair of intersection points and element id
   typedef std::pair<SK::Circular_arc_point_3, UInt> pair_type;
 
 public:
   /// Construct from mesh
   explicit MeshSphereIntersector(Mesh & mesh);
 
   /// Destructor
   virtual ~MeshSphereIntersector();
 
-  /* ------------------------------------------------------------------------ */
-  /* Accessors                                                                */
-  /* ------------------------------------------------------------------------ */
-public:
-  /// get the new_node_per_elem array
-  AKANTU_GET_MACRO(NewNodePerElem, new_node_per_elem, const Array<UInt>)
 
 public:
   /// Construct the primitive tree object
   virtual void constructData();
 
   /**
    * @brief Computes the intersection of the mesh with a sphere
    *
    * @param query (sphere) to compute the intersections with the mesh
    */
   virtual void computeIntersectionQuery(const SK::Sphere_3 & query);
 
+  /// Compute intersection points between the mesh and a query
+  virtual void computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query);
+
   /// Build the IGFEM mesh
   virtual void buildResultFromQueryList(const std::list<SK::Sphere_3> & query);
 
   /// Remove the additionnal nodes
   void removeAdditionnalNodes();
 
+  /// Set the tolerance
+  void setToleranceIntersectionOnNode(UInt tol) {
+    this->tol_intersection_on_node = tol;
+  }
+
 protected:
-  /// new node per element (column 0: number of new nodes, then odd is the intersection node number and even the ID of the sintersected segment)
-  Array<UInt> new_node_per_elem;
+  /// tolerance for which the intersection is considered on the mesh node (relative to the segment lenght)
+  Real tol_intersection_on_node;
 
   /// number of fem nodes in the initial mesh
   const UInt nb_nodes_fem;
 
   /// number of primitive in an element of the template type
   UInt nb_prim_by_el;
 
 };
  
 __END_AKANTU__
 
 #include "mesh_sphere_intersector_tmpl.hh"
 
 #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_HH__
diff --git a/src/geometry/mesh_sphere_intersector_tmpl.hh b/src/geometry/mesh_sphere_intersector_tmpl.hh
index 1110e4e56..462a31ea7 100644
--- a/src/geometry/mesh_sphere_intersector_tmpl.hh
+++ b/src/geometry/mesh_sphere_intersector_tmpl.hh
@@ -1,414 +1,518 @@
 /**
  * @file mesh_sphere_intersector_tmpl.hh
  *
  * @author Clément Roux-Langlois <clement.roux@epfl.ch>
  *
  * @date creation: Wed june 10 2015
  * @date last modification: Wed June 17 2015
  *
  * @brief Computation of mesh intersection with spheres
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2015 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__
 #define __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__
 
 #include "aka_common.hh"
 #include "mesh_geom_common.hh"
 #include "tree_type_helper.hh"
 #include "mesh_sphere_intersector.hh"
 
 __BEGIN_AKANTU__
 
 template<UInt dim, ElementType type>
 MeshSphereIntersector<dim, type>::MeshSphereIntersector(Mesh & mesh):
-  parent_type(mesh),
-  new_node_per_elem(0, 1 + 4*(dim-1)),
+  MeshGeomIntersector<dim, type, Line_arc<SK>, SK::Sphere_3, SK>(mesh),
+  tol_intersection_on_node(1e-10),
   nb_nodes_fem(mesh.getNodes().getSize()),
   nb_prim_by_el(0)
 {
+  this->new_node_per_elem = new Array<UInt>(0, 1 + 4*(dim-1));
+  this->new_nodes = new Array<Real>(0,dim);
   this->constructData();
 
   for(Mesh::type_iterator it = mesh.firstType(dim); it != mesh.lastType(dim); ++it){
 #if defined(AKANTU_IGFEM)
     if(*it == _triangle_3){
       // Addition of the igfem types in the mesh
       this->mesh.addConnectivityType(_igfem_triangle_4, _not_ghost);
       this->mesh.addConnectivityType(_igfem_triangle_4, _ghost);
       this->mesh.addConnectivityType(_igfem_triangle_5, _not_ghost);
       this->mesh.addConnectivityType(_igfem_triangle_5, _ghost);
     }
 
     else if( (*it != _triangle_3) && (*it != _igfem_triangle_4) && (*it != _igfem_triangle_5) ) {
       AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it);
     }
 
     if( (type == _triangle_3) || (type == _igfem_triangle_4) || (type == _igfem_triangle_5) ){
       this->nb_prim_by_el = 3;
+      const_cast<UInt &>(this->nb_seg_by_el) = 3;
     } else {
       AKANTU_DEBUG_ERROR("Not ready for mesh type " << type);
     }
 #else
     if( (*it != _triangle_3) )
       AKANTU_DEBUG_ERROR("Not ready for mesh type " << *it);
 #endif
   }
 }
 
   template<UInt dim, ElementType type>
 MeshSphereIntersector<dim, type>::~MeshSphereIntersector()
 {}
 
 template<UInt dim, ElementType type>
 void MeshSphereIntersector<dim, type>::constructData() {
 
-  this->new_node_per_elem.resize(this->mesh.getConnectivity(type).getSize());
-  this->new_node_per_elem.clear();
+  this->new_node_per_elem->resize(this->mesh.getConnectivity(type).getSize());
+  this->new_node_per_elem->clear();
 
   MeshGeomIntersector<dim, type, Line_arc<SK>, SK::Sphere_3, SK>::constructData();
 }
 
 template<UInt dim, ElementType type>
 void MeshSphereIntersector<dim, type>::computeIntersectionQuery(const SK::Sphere_3 & query) {
   AKANTU_DEBUG_IN();
 
   Array<Real> & nodes = this->mesh.getNodes();
   UInt nb_node = nodes.getSize() ;
+  Array<UInt> & new_node_per_elem = *(this->new_node_per_elem);
 
   // Tolerance for proximity checks should be defined by user
-  Math::setTolerance(1e-10);
+  Math::setTolerance(tol_intersection_on_node);
   typedef boost::variant<pair_type> sk_inter_res;
 
   TreeTypeHelper<Line_arc<Spherical>, Spherical>::const_iterator
     it = this->factory.getPrimitiveList().begin(),
     end= this->factory.getPrimitiveList().end();
 
   NewNodesEvent new_nodes;
   for (; it != end ; ++it) {
     std::list<sk_inter_res> s_results;
     CGAL::intersection(*it, query, std::back_inserter(s_results));
 
     if (s_results.size() == 1) { // just one point
       if (pair_type * pair = boost::get<pair_type>(&s_results.front())) {
         if (pair->second == 1) { // not a point tangent to the sphere
           // Addition of the new node
           Vector<Real> new_node(dim, 0.0);
           Cartesian::Point_3 point(CGAL::to_double(pair->first.x()),
                                    CGAL::to_double(pair->first.y()),
                                    CGAL::to_double(pair->first.z()));
 
           for (UInt i = 0 ; i < dim ; i++) {
             new_node(i) = point[i];
           }
 
           bool is_on_mesh = false, is_new = true;
           // check if we already compute this intersection for a neighboor element
 	  UInt n = nb_nodes_fem-1;
           for (; n < nb_node ; ++n) {
             Array<Real>::vector_iterator existing_node = nodes.begin(dim) + n;
             if (Math::are_vector_equal(dim, new_node.storage(), existing_node->storage())) {
               is_new = false;
               break;
             }
           }
 
             Cartesian::Point_3 source_cgal(CGAL::to_double(it->source().x()),
                                     CGAL::to_double(it->source().y()),
                                     CGAL::to_double(it->source().z()));
             Cartesian::Point_3 target_cgal(CGAL::to_double(it->target().x()),
                                     CGAL::to_double(it->target().y()),
                                     CGAL::to_double(it->target().z()));
 
             Vector<Real> source(dim), target(dim);
             for (UInt i = 0 ; i < dim ; i++) {
               source(i) = source_cgal[i];
               target(i) = target_cgal[i];
             }
 
             // Check if we are close from a node of the segment
             if (Math::are_vector_equal(dim, source.storage(), new_node.storage()) ||
                 Math::are_vector_equal(dim, target.storage(), new_node.storage())) {
               is_on_mesh = true;
               is_new = false;
             }
 
             if (is_new) {
               nodes.push_back(new_node);
               new_nodes.getList().push_back(nb_node);
               nb_node++;
             }
 
             if (!is_on_mesh) {
               new_node_per_elem(it->id(), 0) += 1;
               new_node_per_elem(it->id(), (2 * new_node_per_elem(it->id(), 0)) - 1) = n;
               new_node_per_elem(it->id(), 2 * new_node_per_elem(it->id(), 0)) = it->segId();
             }
 
             else {
               // if intersection is at a node, write node number (in el) in pennultimate position
               if (Math::are_vector_equal(dim, source.storage(), new_node.storage())) {
                 new_node_per_elem(it->id(), (new_node_per_elem.getNbComponent() - 2)) = it->segId() - 1;
               } else {
                 new_node_per_elem(it->id(), (new_node_per_elem.getNbComponent() - 2)) =
                   it->segId() % this->nb_prim_by_el;
               }
             }
           }
         }
       }
   }
 
   if(new_nodes.getList().getSize())
     this->mesh.sendEvent(new_nodes);
 
   AKANTU_DEBUG_OUT();
 }
 
+template<UInt dim, ElementType type>
+void MeshSphereIntersector<dim, type>:: computeMeshQueryIntersectionPoint(const SK::Sphere_3 & query) {
+  /// function to replace computeIntersectionQuery in a more generic geometry module version
+  // The newNodeEvent is not send from this method who only compute the intersection points
+  AKANTU_DEBUG_IN();
+
+  Array<Real> & nodes = this->mesh.getNodes();
+  UInt nb_node = nodes.getSize() + this->new_nodes->getSize();
+  Array<UInt> & new_node_per_elem = *(this->new_node_per_elem);
+
+  // Tolerance for proximity checks should be defined by user
+  Math::setTolerance(tol_intersection_on_node);
+  typedef boost::variant<pair_type> sk_inter_res;
+
+  TreeTypeHelper<Line_arc<Spherical>, Spherical>::const_iterator
+    it = this->factory.getPrimitiveList().begin(),
+    end= this->factory.getPrimitiveList().end();
+
+  //NewNodesEvent new_nodes;
+  for (; it != end ; ++it) {
+    std::list<sk_inter_res> s_results;
+    CGAL::intersection(*it, query, std::back_inserter(s_results));
+
+    if (s_results.size() == 1) { // just one point
+      if (pair_type * pair = boost::get<pair_type>(&s_results.front())) {
+        if (pair->second == 1) { // not a point tangent to the sphere
+          // Addition of the new node
+          Vector<Real> new_node(dim, 0.0);
+          Cartesian::Point_3 point(CGAL::to_double(pair->first.x()),
+                                   CGAL::to_double(pair->first.y()),
+                                   CGAL::to_double(pair->first.z()));
+
+          for (UInt i = 0 ; i < dim ; i++) {
+            new_node(i) = point[i];
+          }
+
+          bool is_on_mesh = false, is_new = true;
+          // check if we already compute this intersection for a neighboor element
+	  UInt n = nb_nodes_fem-1;
+          for (; n < nb_node ; ++n) {
+            Array<Real>::vector_iterator existing_node = nodes.begin(dim) + n;
+            if (Math::are_vector_equal(dim, new_node.storage(), existing_node->storage())) {
+              is_new = false;
+              break;
+            }
+          }
+
+            Cartesian::Point_3 source_cgal(CGAL::to_double(it->source().x()),
+                                    CGAL::to_double(it->source().y()),
+                                    CGAL::to_double(it->source().z()));
+            Cartesian::Point_3 target_cgal(CGAL::to_double(it->target().x()),
+                                    CGAL::to_double(it->target().y()),
+                                    CGAL::to_double(it->target().z()));
+
+            Vector<Real> source(dim), target(dim);
+            for (UInt i = 0 ; i < dim ; i++) {
+              source(i) = source_cgal[i];
+              target(i) = target_cgal[i];
+            }
+
+            // Check if we are close from a node of the segment
+            if (Math::are_vector_equal(dim, source.storage(), new_node.storage()) ||
+                Math::are_vector_equal(dim, target.storage(), new_node.storage())) {
+              is_on_mesh = true;
+              is_new = false;
+            }
+
+            if (is_new) {
+              this->new_nodes->push_back(new_node);
+              // new_nodes.getList().push_back(nb_node);
+              nb_node++;
+            }
+
+            if (!is_on_mesh) {
+              new_node_per_elem(it->id(), 0) += 1;
+              new_node_per_elem(it->id(), (2 * new_node_per_elem(it->id(), 0)) - 1) = n;
+              new_node_per_elem(it->id(), 2 * new_node_per_elem(it->id(), 0)) = it->segId();
+            }
+
+            else {
+              // if intersection is at a node, write node number (in el) in pennultimate position
+              if (Math::are_vector_equal(dim, source.storage(), new_node.storage())) {
+                new_node_per_elem(it->id(), (new_node_per_elem.getNbComponent() - 2)) = it->segId() - 1;
+              } else {
+                new_node_per_elem(it->id(), (new_node_per_elem.getNbComponent() - 2)) =
+                  it->segId() % this->nb_prim_by_el;
+              }
+            }
+          }
+        }
+      }
+  }
+
+  /*if(new_nodes.getList().getSize())
+    this->mesh.sendEvent(new_nodes);*/
+
+  AKANTU_DEBUG_OUT();
+}
+
 template<UInt dim, ElementType type>
 void MeshSphereIntersector<dim, type>::removeAdditionnalNodes() {
   AKANTU_DEBUG_IN();
 
   RemovedNodesEvent remove_nodes(this->mesh);
   Array<UInt> & nodes_removed = remove_nodes.getList();
   Array<UInt> & new_numbering = remove_nodes.getNewNumbering();
 
   for(UInt i = 0 ; i < nb_nodes_fem ; ++i){
     new_numbering(i) = i;
   }
 
   for(UInt i = nb_nodes_fem ; i < this->mesh.getNodes().getSize(); ++i){
     new_numbering(i) = UInt(-1) ;
     nodes_removed.push_back(i);
   }
 
   this->mesh.sendEvent(remove_nodes);
   AKANTU_DEBUG_OUT();
 }
 
 
 #if defined(AKANTU_IGFEM)
 
 template<UInt dim, ElementType type>
 void MeshSphereIntersector<dim, type>::buildResultFromQueryList(const std::list<SK::Sphere_3> & query_list) {
   AKANTU_DEBUG_IN();
 
   this->computeIntersectionQueryList(query_list);
 
   Array<UInt> & connec_type_tmpl = this->mesh.getConnectivity(type);
+  const Array<UInt> new_node_per_elem = *(this->new_node_per_elem);
 
   Array<UInt>
     & connec_igfem_tri4 = this->mesh.getConnectivity(_igfem_triangle_4),
     & connec_igfem_tri5 = this->mesh.getConnectivity(_igfem_triangle_5),
     & connec_tri3 = this->mesh.getConnectivity(_triangle_3);
 
   Element element_tri3(_triangle_3, 0, _not_ghost, _ek_regular);
   Element element_tri4(_igfem_triangle_4, 0, _not_ghost, _ek_igfem);
   Element element_tri5(_igfem_triangle_5, 0, _not_ghost, _ek_igfem);
 
   NewIGFEMElementsEvent new_elements;
   UInt n_new_el = 0;
   Array<Element> & new_elements_list = new_elements.getList();
   Array<Element> & old_elements_list = new_elements.getOldElementsList();
  
   RemovedElementsEvent remove_elem(this->mesh);
   remove_elem.getNewNumbering().alloc(connec_type_tmpl.getSize(), 1, type, _not_ghost);
   Array<UInt> & new_numbering = remove_elem.getNewNumbering(type, _not_ghost);
   UInt new_nb_type_tmpl = 0; // Meter of the number of element (type) will we loop on elements
 
   for (UInt nel = 0 ; nel != new_node_per_elem.getSize(); ++nel) {
     if( (type != _triangle_3)
         && ( (new_node_per_elem(nel,0)==0)
              || ( (new_node_per_elem(nel,0) == 1) 
                   && ( ( (new_node_per_elem(nel,2)+1) % this->nb_prim_by_el ) 
                        != new_node_per_elem(nel, new_node_per_elem.getNbComponent() - 2) ) ) ) ){
       Element element_type_tmpl(type, 0, _not_ghost, Mesh::getKind(type));  
       new_elements_list.resize(n_new_el+1);
       old_elements_list.resize(n_new_el+1);
       Vector<UInt> ctri3(3);
       ctri3(0) =  connec_type_tmpl(nel,0);
       ctri3(1) =  connec_type_tmpl(nel,1);
       ctri3(2) =  connec_type_tmpl(nel,2);
       /// add the new element in the mesh
       UInt nb_tri3 = connec_tri3.getSize();
       connec_tri3.push_back(ctri3);
       element_tri3.element = nb_tri3;
       new_elements_list(n_new_el) = element_tri3;
       /// the number of the old element in the mesh
       element_type_tmpl.element = nel;
       old_elements_list(n_new_el) = element_type_tmpl;
       new_numbering(nel) =  UInt(-1);
       ++n_new_el;
 
       remove_elem.getList().push_back(element_type_tmpl);
   
     }
     else if( (new_node_per_elem(nel,0)!=0)
              && !( (new_node_per_elem(nel,0) == 1) 
                    && ( ( (new_node_per_elem(nel,2)+1) % this->nb_prim_by_el ) 
                         != new_node_per_elem(nel, new_node_per_elem.getNbComponent() - 2) ) ) ){
       Element element_type_tmpl(type, 0, _not_ghost);
       element_type_tmpl.kind = Mesh::getKind(type);
       new_elements_list.resize(n_new_el+1);
       old_elements_list.resize(n_new_el+1);
       switch(new_node_per_elem(nel,0)){
       case 1 :{
 	Vector<UInt> ctri4(4);
 	switch(new_node_per_elem(nel,2)){
 	case 1 :
 	  ctri4(0) = connec_type_tmpl(nel,2);
 	  ctri4(1) = connec_type_tmpl(nel,0);
 	  ctri4(2) = connec_type_tmpl(nel,1);
 	  break;
 	case 2 :
 	  ctri4(0) = connec_type_tmpl(nel,0);
 	  ctri4(1) = connec_type_tmpl(nel,1);
 	  ctri4(2) = connec_type_tmpl(nel,2);
 	  break;
 	case 3 :
 	  ctri4(0) = connec_type_tmpl(nel,1);
 	  ctri4(1) = connec_type_tmpl(nel,2);
 	  ctri4(2) = connec_type_tmpl(nel,0);
 	  break;
 	default :
 	  AKANTU_DEBUG_ERROR("A triangle have only 3 segments not "<< new_node_per_elem(nel,2));
 	  break;
 	}
 	ctri4(3) = new_node_per_elem(nel,1);
 	UInt nb_tri4 = connec_igfem_tri4.getSize();
 	connec_igfem_tri4.push_back(ctri4);
 	element_tri4.element = nb_tri4;
 	//new_elements.getList().push_back(element_tri4);
 	new_elements_list(n_new_el) = element_tri4;
 	if(type == _igfem_triangle_4)
 	   new_numbering.push_back(connec_igfem_tri4.getSize()-2);
 	break;
       }
       case 2 :{
 	Vector<UInt> ctri5(5);
 	if( (new_node_per_elem(nel,2)==1) && (new_node_per_elem(nel,4)==2) ){
 	  ctri5(0) = connec_type_tmpl(nel,1);
 	  ctri5(1) = connec_type_tmpl(nel,2);
 	  ctri5(2) = connec_type_tmpl(nel,0);
 	  ctri5(3) = new_node_per_elem(nel,3);
 	  ctri5(4) = new_node_per_elem(nel,1);
 	}
 	else if((new_node_per_elem(nel,2)==1) && (new_node_per_elem(nel,4)==3)){
 	  ctri5(0) = connec_type_tmpl(nel,0);
 	  ctri5(1) = connec_type_tmpl(nel,1);
 	  ctri5(2) = connec_type_tmpl(nel,2);
 	  ctri5(3) = new_node_per_elem(nel,1);
 	  ctri5(4) = new_node_per_elem(nel,3);
 	}
 	else if((new_node_per_elem(nel,2)==2) && (new_node_per_elem(nel,4)==3)){
 	  ctri5(0) = connec_type_tmpl(nel,2);
 	  ctri5(1) = connec_type_tmpl(nel,0);
 	  ctri5(2) = connec_type_tmpl(nel,1);
 	  ctri5(3) = new_node_per_elem(nel,3);
 	  ctri5(4) = new_node_per_elem(nel,1);
 	}
 	else if((new_node_per_elem(nel,2)==2) && (new_node_per_elem(nel,4)==1)){
 	  ctri5(0) = connec_type_tmpl(nel,1);
 	  ctri5(1) = connec_type_tmpl(nel,2);
 	  ctri5(2) = connec_type_tmpl(nel,0);
 	  ctri5(3) = new_node_per_elem(nel,1);
 	  ctri5(4) = new_node_per_elem(nel,3);
 	}
 	else if((new_node_per_elem(nel,2)==3) && (new_node_per_elem(nel,4)==1)){
 	  ctri5(0) = connec_type_tmpl(nel,0);
 	  ctri5(1) = connec_type_tmpl(nel,1);
 	  ctri5(2) = connec_type_tmpl(nel,2);
 	  ctri5(3) = new_node_per_elem(nel,3);
 	  ctri5(4) = new_node_per_elem(nel,1);
 	}
 	else if((new_node_per_elem(nel,2)==3) && (new_node_per_elem(nel,4)==2)){
 	  ctri5(0) = connec_type_tmpl(nel,2);
 	  ctri5(1) = connec_type_tmpl(nel,0);
 	  ctri5(2) = connec_type_tmpl(nel,1);
 	  ctri5(3) = new_node_per_elem(nel,1);
 	  ctri5(4) = new_node_per_elem(nel,3);
 	}
 	else{
 	  AKANTU_DEBUG_ERROR("A triangle have only 3 segments not "<< new_node_per_elem(nel,2) << "and" << new_node_per_elem(nel,4));
 	}	      
 	UInt nb_tri5 = connec_igfem_tri5.getSize();
 	connec_igfem_tri5.push_back(ctri5);
 	element_tri5.element = nb_tri5;
 	//new_elements.getList().push_back(element_tri5);
 	new_elements_list(n_new_el) = element_tri5;
 	if(type == _igfem_triangle_5){
 	  new_numbering.push_back(new_nb_type_tmpl);
 	  ++new_nb_type_tmpl;
 	}
 	break;
       }
       default:
 	AKANTU_DEBUG_ERROR("Igfem cannot add "<< new_node_per_elem(nel,0) << " nodes to triangles");
 	break;
       }
       element_type_tmpl.element = nel;
       old_elements_list(n_new_el) = element_type_tmpl;
       remove_elem.getList().push_back(element_type_tmpl);
       new_numbering(nel) =  UInt(-1);
       ++n_new_el;
     }
     else{
       new_numbering(nel) = new_nb_type_tmpl;
       ++new_nb_type_tmpl;
     }
   }
 
   // for(UInt nel=new_node_per_elem.getSize(); nel < this->mesh.getNbElement(type); ++nel) {
   //     new_numbering(nel) = new_nb_type_tmpl;
   //     ++new_nb_type_tmpl;
   // }
 
   UInt el_index = 0;
   for (UInt e = 0; e < this->mesh.getNbElement(type); ++e) {
     if (new_numbering(e) != UInt(-1)) {
       new_numbering(e) = el_index;
       ++el_index;
     }
   }
     
   if(n_new_el > 0){
     this->mesh.sendEvent(new_elements);
     this->mesh.sendEvent(remove_elem);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 #else
 
 template<UInt dim, ElementType type>
 void MeshSphereIntersector<dim, type>::buildResultFromQueryList(const std::list<SK::Sphere_3> & query_list) {
   AKANTU_DEBUG_TO_IMPLEMENT();
 }
 
 
 #endif
 
 __END_AKANTU__
 
 #endif // __AKANTU_MESH_SPHERE_INTERSECTOR_TMPL_HH__
 
diff --git a/third-party/cmake/blackdynamite.cmake b/third-party/cmake/blackdynamite.cmake
index 868a3c6e6..e37bcdaa6 100644
--- a/third-party/cmake/blackdynamite.cmake
+++ b/third-party/cmake/blackdynamite.cmake
@@ -1,41 +1,41 @@
 if(${PROJECT_SOURCE_DIR}/third-party/${BLACKDYNAMITE_ARCHIVE})
   set(_blackdynamite_download_command
     URL ${PROJECT_SOURCE_DIR}/third-party/${BLACKDYNAMITE_ARCHIVE})
 else()
   set(_blackdynamite_download_command
-    SVN_REPOSITORY ${BLACKDYNAMITE_URL}
+    GIT_REPOSITORY ${BLACKDYNAMITE_GIT}
+    GIT_TAG ${BLACKDYNAMITE_VERSION}
     )
 endif()
 
 if(CMAKE_VERSION VERSION_GREATER 3.1)
   set(_extra_options 
     UPDATE_DISCONNECTED 1
     DOWNLOAD_NO_PROGRESS 1
     EXCLUDE_FROM_ALL 1
     )
 endif()
 
 ExternalProject_Add(blackdynamite
   PREFIX ${PROJECT_BINARY_DIR}/third-party
   ${_blackdynamite_download_command}
   ${_extra_options}
   CMAKE_ARGS <SOURCE_DIR>/
   CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=<INSTALL_DIR> -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE} -DCMAKE_CXX_COMPILER:PATH=${CMAKE_CXX_COMPILER}
   BUILD_COMMAND make
   INSTALL_COMMAND make install
   LOG_DOWNLOAD 1
   LOG_CONFIGURE 1
   LOG_BUILD 1
   LOG_INSTALL 1
   )
 
 
-
 set_third_party_shared_libirary_name(BLACKDYNAMITE_LIBRARIES blackdynamite)
-package_set_include_dir(BLACKDYNAMITE_INCLUDE_DIR "${PROJECT_BINARY_DIR}/third-party/include/blackdynamite" CACHE PATH "")
+set(BLACKDYNAMITE_INCLUDE_DIR "${PROJECT_BINARY_DIR}/third-party/include/blackdynamite" CACHE PATH "")
 mark_as_advanced(
   BLACKDYNAMITE_LIBRARIES
   BLACKDYNAMITE_INCLUDE_DIR
   )
 
 package_add_extra_dependency(BlackDynamite blackdynamite)