diff --git a/packages/core.cmake b/packages/core.cmake
index 10660f10e..c8d23637e 100644
--- a/packages/core.cmake
+++ b/packages/core.cmake
@@ -1,433 +1,440 @@
 #===============================================================================
 # @file   core.cmake
 #
 # @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @date creation: Mon Nov 21 2011
 # @date last modification: Fri Apr 09 2021
 #
 # @brief  package description for core
 #
 #
 # @section LICENSE
 #
 # Copyright (©) 2010-2021 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(core NOT_OPTIONAL
   DESCRIPTION "core package for Akantu"
   FEATURES_PUBLIC cxx_strong_enums cxx_defaulted_functions
                   cxx_deleted_functions cxx_auto_type cxx_decltype_auto
   FEATURES_PRIVATE cxx_lambdas cxx_nullptr cxx_range_for
                    cxx_delegating_constructors
   DEPENDS INTERFACE akantu_iterators Boost
   )
 
 if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU" OR CMAKE_CXX_COMPILER_ID STREQUAL "Clang")
   package_set_compile_flags(core "-Wall -Wextra -pedantic")
 else()
   package_set_compile_flags(core "-Wall")
 endif()
 
 package_declare_sources(core
   common/aka_array.cc
   common/aka_array.hh
   common/aka_array_filter.hh
   common/aka_array_tmpl.hh
   common/aka_array_printer.hh
   common/aka_bbox.hh
   common/aka_blas_lapack.hh
   common/aka_circular_array.hh
   common/aka_circular_array_inline_impl.hh
   common/aka_common.cc
   common/aka_common.hh
   common/aka_common_inline_impl.hh
   common/aka_csr.hh
   common/aka_element_classes_info_inline_impl.hh
   common/aka_enum_macros.hh
   common/aka_error.cc
   common/aka_error.hh
   common/aka_event_handler_manager.hh
   common/aka_extern.cc
   common/aka_factory.hh
   common/aka_fwd.hh
   common/aka_grid_dynamic.hh
   common/aka_math.cc
   common/aka_math.hh
   common/aka_math_tmpl.hh
   common/aka_named_argument.hh
   common/aka_random_generator.hh
   common/aka_safe_enum.hh
   common/aka_types.hh
   common/aka_voigthelper.hh
   common/aka_voigthelper_tmpl.hh
   common/aka_voigthelper.cc
   common/aka_warning.hh
   common/aka_warning_restore.hh
   
   fe_engine/element_class.hh
   fe_engine/element_class_helper.hh
   fe_engine/element_class_tmpl.hh
   fe_engine/element_classes/element_class_hexahedron_8_inline_impl.hh
   fe_engine/element_classes/element_class_hexahedron_20_inline_impl.hh
   fe_engine/element_classes/element_class_pentahedron_6_inline_impl.hh
   fe_engine/element_classes/element_class_pentahedron_15_inline_impl.hh
   fe_engine/element_classes/element_class_point_1_inline_impl.hh
   fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh
   fe_engine/element_classes/element_class_quadrangle_8_inline_impl.hh
   fe_engine/element_classes/element_class_segment_2_inline_impl.hh
   fe_engine/element_classes/element_class_segment_3_inline_impl.hh
   fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.hh
   fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.hh
   fe_engine/element_classes/element_class_triangle_3_inline_impl.hh
   fe_engine/element_classes/element_class_triangle_6_inline_impl.hh
   fe_engine/element_type_conversion.hh
 
   fe_engine/fe_engine.cc
   fe_engine/fe_engine.hh
   fe_engine/fe_engine_inline_impl.hh
   fe_engine/fe_engine_template.hh
   fe_engine/fe_engine_template_tmpl_field.hh
   fe_engine/fe_engine_template_tmpl.hh
   fe_engine/geometrical_element_property.hh
   fe_engine/geometrical_element_property.cc
   fe_engine/gauss_integration.cc
   fe_engine/gauss_integration_tmpl.hh
   fe_engine/integrator.hh
   fe_engine/integrator_gauss.hh
   fe_engine/integrator_gauss_inline_impl.hh
   fe_engine/interpolation_element_tmpl.hh
   fe_engine/integration_point.hh
   fe_engine/shape_functions.hh
   fe_engine/shape_functions.cc
   fe_engine/shape_functions_inline_impl.hh
   fe_engine/shape_lagrange_base.cc
   fe_engine/shape_lagrange_base.hh
   fe_engine/shape_lagrange_base_inline_impl.hh
   fe_engine/shape_lagrange.hh
   fe_engine/shape_lagrange_inline_impl.hh
   fe_engine/element.hh
 
   io/dumper/dumpable.hh
   io/dumper/dumpable.cc
   io/dumper/dumpable_dummy.hh
   io/dumper/dumpable_inline_impl.hh
   io/dumper/dumper_field.hh
   io/dumper/dumper_material_padders.hh
   io/dumper/dumper_filtered_connectivity.hh
   io/dumper/dumper_element_partition.hh
 
   io/mesh_io.cc
   io/mesh_io.hh
   io/mesh_io/mesh_io_diana.cc
   io/mesh_io/mesh_io_diana.hh
   io/mesh_io/mesh_io_msh.cc
   io/mesh_io/mesh_io_msh.hh
-  #io/model_io.cc
-  #io/model_io.hh
 
   io/parser/algebraic_parser.hh
   io/parser/input_file_parser.hh
   io/parser/parsable.cc
   io/parser/parsable.hh
   io/parser/parser.cc
   io/parser/parser_real.cc
   io/parser/parser_random.cc
   io/parser/parser_types.cc
   io/parser/parser_input_files.cc
   io/parser/parser.hh
   io/parser/parser_tmpl.hh
   io/parser/parser_grammar_tmpl.hh
   io/parser/cppargparse/cppargparse.hh
   io/parser/cppargparse/cppargparse.cc
   io/parser/cppargparse/cppargparse_tmpl.hh
 
   io/parser/parameter_registry.cc
   io/parser/parameter_registry.hh
   io/parser/parameter_registry_tmpl.hh
 
   mesh/element_group.cc
   mesh/element_group.hh
   mesh/element_group_inline_impl.hh
   mesh/element_type_map.cc
   mesh/element_type_map.hh
   mesh/element_type_map_tmpl.hh
   mesh/element_type_map_filter.hh
   mesh/group_manager.cc
   mesh/group_manager.hh
   mesh/group_manager_inline_impl.hh
   mesh/mesh.cc
   mesh/mesh.hh
   mesh/mesh_periodic.cc
   mesh/mesh_accessor.hh
   mesh/mesh_events.hh
   mesh/mesh_filter.hh
   mesh/mesh_global_data_updater.hh
   mesh/mesh_data.cc
   mesh/mesh_data.hh
   mesh/mesh_data_tmpl.hh
   mesh/mesh_inline_impl.hh
   mesh/node_group.cc
   mesh/node_group.hh
   mesh/node_group_inline_impl.hh
   mesh/mesh_iterators.hh
 
   mesh_utils/mesh_partition.cc
   mesh_utils/mesh_partition.hh
   mesh_utils/mesh_partition/mesh_partition_mesh_data.cc
   mesh_utils/mesh_partition/mesh_partition_mesh_data.hh
   mesh_utils/mesh_partition/mesh_partition_scotch.hh
   mesh_utils/mesh_utils_pbc.cc
   mesh_utils/mesh_utils.cc
   mesh_utils/mesh_utils.hh
   mesh_utils/mesh_utils_distribution.cc
   mesh_utils/mesh_utils_distribution.hh
   mesh_utils/mesh_utils.hh
   mesh_utils/mesh_utils_inline_impl.hh
   mesh_utils/global_ids_updater.hh
   mesh_utils/global_ids_updater.cc
   mesh_utils/global_ids_updater_inline_impl.hh
 
   model/common/boundary_condition/boundary_condition.hh
   model/common/boundary_condition/boundary_condition_functor.hh
   model/common/boundary_condition/boundary_condition_functor_inline_impl.hh
   model/common/boundary_condition/boundary_condition_tmpl.hh
 
+  model/common/constitutive_laws/constitutive_law_non_local_interface.hh
+  model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
+  model/common/constitutive_laws/constitutive_laws_handler.hh
+  model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh
+  model/common/constitutive_laws/internal_field.hh
+  model/common/constitutive_laws/internal_field_tmpl.hh
+  model/common/constitutive_laws/random_internal_field.hh
+  model/common/constitutive_laws/random_internal_field_tmpl.hh
+
   model/common/non_local_toolbox/neighborhood_base.hh
   model/common/non_local_toolbox/neighborhood_base.cc
   model/common/non_local_toolbox/neighborhood_base_inline_impl.hh
   model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh
   model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc
   model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.hh
   model/common/non_local_toolbox/non_local_manager.hh
   model/common/non_local_toolbox/non_local_manager.cc
   model/common/non_local_toolbox/non_local_manager_inline_impl.hh
   model/common/non_local_toolbox/non_local_manager_callback.hh
   model/common/non_local_toolbox/non_local_neighborhood_base.hh
   model/common/non_local_toolbox/non_local_neighborhood_base.cc
   model/common/non_local_toolbox/non_local_neighborhood.hh
   model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh
   model/common/non_local_toolbox/non_local_neighborhood_inline_impl.hh
   model/common/non_local_toolbox/base_weight_function.hh
   model/common/non_local_toolbox/base_weight_function_inline_impl.hh
 
   model/common/model_solver.cc
   model/common/model_solver.hh
   model/common/solver_callback.hh
   model/common/solver_callback.cc
 
   model/common/dof_manager/dof_manager.cc
   model/common/dof_manager/dof_manager.hh
   model/common/dof_manager/dof_manager_default.cc
   model/common/dof_manager/dof_manager_default.hh
   model/common/dof_manager/dof_manager_default_inline_impl.hh
   model/common/dof_manager/dof_manager_inline_impl.hh
 
   model/common/non_linear_solver/non_linear_solver.cc
   model/common/non_linear_solver/non_linear_solver.hh
   model/common/non_linear_solver/non_linear_solver_default.hh
   model/common/non_linear_solver/non_linear_solver_lumped.cc
   model/common/non_linear_solver/non_linear_solver_lumped.hh
 
   model/common/time_step_solvers/time_step_solver.hh
   model/common/time_step_solvers/time_step_solver.cc
   model/common/time_step_solvers/time_step_solver_default.cc
   model/common/time_step_solvers/time_step_solver_default.hh
   model/common/time_step_solvers/time_step_solver_default_explicit.hh
 
   model/common/integration_scheme/generalized_trapezoidal.cc
   model/common/integration_scheme/generalized_trapezoidal.hh
   model/common/integration_scheme/integration_scheme.cc
   model/common/integration_scheme/integration_scheme.hh
   model/common/integration_scheme/integration_scheme_1st_order.cc
   model/common/integration_scheme/integration_scheme_1st_order.hh
   model/common/integration_scheme/integration_scheme_2nd_order.cc
   model/common/integration_scheme/integration_scheme_2nd_order.hh
   model/common/integration_scheme/newmark-beta.cc
   model/common/integration_scheme/newmark-beta.hh
   model/common/integration_scheme/pseudo_time.cc
   model/common/integration_scheme/pseudo_time.hh
 
   model/model.cc
   model/model.hh
   model/model_inline_impl.hh
   model/model_options.hh
 
   solver/solver_vector.hh
   solver/solver_vector_default.hh
   solver/solver_vector_default_tmpl.hh
   solver/solver_vector_distributed.cc
   solver/solver_vector_distributed.hh
   solver/sparse_matrix.cc
   solver/sparse_matrix.hh
   solver/sparse_matrix_aij.cc
   solver/sparse_matrix_aij.hh
   solver/sparse_matrix_aij_inline_impl.hh
   solver/sparse_matrix_inline_impl.hh
   solver/sparse_solver.cc
   solver/sparse_solver.hh
   solver/sparse_solver_inline_impl.hh
   solver/terms_to_assemble.hh
 
   synchronizer/communication_buffer_inline_impl.hh
   synchronizer/communication_descriptor.hh
   synchronizer/communication_descriptor_tmpl.hh
   synchronizer/communication_request.hh
   synchronizer/communication_tag.hh
   synchronizer/communications.hh
   synchronizer/communications_tmpl.hh
   synchronizer/communicator.cc
   synchronizer/communicator.hh
   synchronizer/communicator_dummy_inline_impl.hh
   synchronizer/communicator_event_handler.hh
   synchronizer/communicator_inline_impl.hh
   synchronizer/data_accessor.cc
   synchronizer/data_accessor.hh
   synchronizer/dof_synchronizer.cc
   synchronizer/dof_synchronizer.hh
   synchronizer/dof_synchronizer_inline_impl.hh
   synchronizer/element_info_per_processor.cc
   synchronizer/element_info_per_processor.hh
   synchronizer/element_info_per_processor_tmpl.hh
   synchronizer/element_synchronizer.cc
   synchronizer/element_synchronizer.hh
   synchronizer/facet_synchronizer.cc
   synchronizer/facet_synchronizer.hh
   synchronizer/facet_synchronizer_inline_impl.hh
   synchronizer/grid_synchronizer.cc
   synchronizer/grid_synchronizer.hh
   synchronizer/grid_synchronizer_tmpl.hh
   synchronizer/master_element_info_per_processor.cc
   synchronizer/node_info_per_processor.cc
   synchronizer/node_info_per_processor.hh
   synchronizer/node_synchronizer.cc
   synchronizer/node_synchronizer.hh
   synchronizer/node_synchronizer_inline_impl.hh
   synchronizer/periodic_node_synchronizer.cc
   synchronizer/periodic_node_synchronizer.hh
   synchronizer/slave_element_info_per_processor.cc
   synchronizer/synchronizer.cc
   synchronizer/synchronizer.hh
   synchronizer/synchronizer_impl.hh
   synchronizer/synchronizer_impl_tmpl.hh
   synchronizer/synchronizer_registry.cc
   synchronizer/synchronizer_registry.hh
   synchronizer/synchronizer_tmpl.hh
   synchronizer/communication_buffer.hh
   )
 
 set(AKANTU_SPIRIT_SOURCES
   io/mesh_io/mesh_io_abaqus.cc
   io/parser/parser_real.cc
   io/parser/parser_random.cc
   io/parser/parser_types.cc
   io/parser/parser_input_files.cc
   PARENT_SCOPE
   )
 
 package_declare_elements(core
   ELEMENT_TYPES
   _point_1
   _segment_2
   _segment_3
   _triangle_3
   _triangle_6
   _quadrangle_4
   _quadrangle_8
   _tetrahedron_4
   _tetrahedron_10
   _pentahedron_6
   _pentahedron_15
   _hexahedron_8
   _hexahedron_20
   KIND regular
   GEOMETRICAL_TYPES
   _gt_point
   _gt_segment_2
   _gt_segment_3
   _gt_triangle_3
   _gt_triangle_6
   _gt_quadrangle_4
   _gt_quadrangle_8
   _gt_tetrahedron_4
   _gt_tetrahedron_10
   _gt_hexahedron_8
   _gt_hexahedron_20
   _gt_pentahedron_6
   _gt_pentahedron_15
   INTERPOLATION_TYPES
   _itp_lagrange_point_1
   _itp_lagrange_segment_2
   _itp_lagrange_segment_3
   _itp_lagrange_triangle_3
   _itp_lagrange_triangle_6
   _itp_lagrange_quadrangle_4
   _itp_serendip_quadrangle_8
   _itp_lagrange_tetrahedron_4
   _itp_lagrange_tetrahedron_10
   _itp_lagrange_hexahedron_8
   _itp_serendip_hexahedron_20
   _itp_lagrange_pentahedron_6
   _itp_lagrange_pentahedron_15
   GEOMETRICAL_SHAPES
   _gst_point
   _gst_triangle
   _gst_square
   _gst_prism
   GAUSS_INTEGRATION_TYPES
   _git_point
   _git_segment
   _git_triangle
   _git_tetrahedron
   _git_pentahedron
   INTERPOLATION_KIND _itk_lagrangian
   FE_ENGINE_LISTS
   gradient_on_integration_points
   interpolate_on_integration_points
   interpolate
   compute_normals_on_integration_points
   inverse_map
   contains
   compute_shapes
   compute_shapes_derivatives
   get_N
   compute_dnds
   compute_d2nds2
   compute_jmat
   get_shapes_derivatives
   lagrange_base
   assemble_fields
   )
 
 find_program(READLINK_COMMAND readlink)
 find_program(ADDR2LINE_COMMAND addr2line)
 find_program(PATCH_COMMAND patch)
 mark_as_advanced(READLINK_COMMAND)
 mark_as_advanced(ADDR2LINE_COMMAND)
 
 package_declare_extra_files_to_package(core
   SOURCES
     common/aka_element_classes_info.hh.in
     common/aka_config.hh.in
   )
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang" AND (NOT CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.9))
   package_set_compile_flags(core CXX "-Wno-undefined-var-template")
 endif()
diff --git a/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake
index c7b8776ac..4b86b201f 100644
--- a/packages/solid_mechanics.cmake
+++ b/packages/solid_mechanics.cmake
@@ -1,118 +1,110 @@
 #===============================================================================
 # @file   solid_mechanics.cmake
 #
 # @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
 # @author Nicolas Richart <nicolas.richart@epfl.ch>
 #
 # @date creation: Mon Dec 04 2017
 # @date last modification: Fri Mar 26 2021
 #
 # @brief  package description for core
 #
 #
 # @section LICENSE
 #
 # Copyright (©) 2016-2021 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(solid_mechanics DEFAULT ON
   DESCRIPTION "Solid mechanics model"
   DEPENDS core lapack
   )
 
 package_declare_sources(solid_mechanics
   model/solid_mechanics/material.cc
   model/solid_mechanics/material.hh
   model/solid_mechanics/material_inline_impl.hh
   model/solid_mechanics/material_selector.hh
   model/solid_mechanics/material_selector_tmpl.hh
-  model/solid_mechanics/materials/internal_field.hh
-  model/solid_mechanics/materials/internal_field_tmpl.hh
-  model/solid_mechanics/materials/random_internal_field.hh
-  model/solid_mechanics/materials/random_internal_field_tmpl.hh
   model/solid_mechanics/solid_mechanics_model.cc
   model/solid_mechanics/solid_mechanics_model.hh
-  model/solid_mechanics/solid_mechanics_model_inline_impl.hh
   model/solid_mechanics/solid_mechanics_model_io.cc
   model/solid_mechanics/solid_mechanics_model_mass.cc
-  model/solid_mechanics/solid_mechanics_model_material.cc
-  model/solid_mechanics/solid_mechanics_model_tmpl.hh
   model/solid_mechanics/solid_mechanics_model_event_handler.hh
   model/solid_mechanics/materials/plane_stress_toolbox.hh
   model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh
 
   model/solid_mechanics/materials/material_core_includes.hh
   model/solid_mechanics/materials/material_elastic.cc
   model/solid_mechanics/materials/material_elastic.hh
   model/solid_mechanics/materials/material_elastic_inline_impl.hh
   model/solid_mechanics/materials/material_thermal.cc
   model/solid_mechanics/materials/material_thermal.hh
   model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
   model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
   model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.hh
   model/solid_mechanics/materials/material_elastic_orthotropic.cc
   model/solid_mechanics/materials/material_elastic_orthotropic.hh
   model/solid_mechanics/materials/material_damage/material_anisotropic_damage.hh
   model/solid_mechanics/materials/material_damage/material_anisotropic_damage.cc
   model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh
   model/solid_mechanics/materials/material_damage/material_damage.hh
   model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh
   model/solid_mechanics/materials/material_damage/material_marigo.cc
   model/solid_mechanics/materials/material_damage/material_marigo.hh
   model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.hh
   model/solid_mechanics/materials/material_damage/material_mazars.cc
   model/solid_mechanics/materials/material_damage/material_mazars.hh
   model/solid_mechanics/materials/material_damage/material_phasefield.cc
   model/solid_mechanics/materials/material_damage/material_phasefield.hh
   model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh
   model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.hh
   model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc
   model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh
   model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.hh
   model/solid_mechanics/materials/material_plastic/material_plastic.cc
   model/solid_mechanics/materials/material_plastic/material_plastic.hh
   model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.hh
   model/solid_mechanics/materials/material_plastic/material_drucker_prager.cc
   model/solid_mechanics/materials/material_plastic/material_drucker_prager.hh
   model/solid_mechanics/materials/material_plastic/material_drucker_prager_inline_impl.hh
   model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc
   model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh
   model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.hh
   model/solid_mechanics/materials/material_damage/material_von_mises_mazars.cc
   model/solid_mechanics/materials/material_damage/material_von_mises_mazars.hh
   model/solid_mechanics/materials/material_damage/material_von_mises_mazars_inline_impl.hh
   model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc
   model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh
   model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc
   model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh
 
   model/solid_mechanics/materials/material_non_local.hh
-  model/solid_mechanics/materials/material_non_local_tmpl.hh
   model/solid_mechanics/materials/material_non_local_includes.hh
   )
 
 package_declare_material_infos(solid_mechanics
   LIST AKANTU_CORE_MATERIAL_LIST
   INCLUDE material_core_includes.hh
   )
 
 package_declare_extra_files_to_package(solid_mechanics
   SOURCES
     model/solid_mechanics/material_list.hh.in
   )
diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh
index fa1966f8f..3302f921f 100644
--- a/src/io/parser/parser.hh
+++ b/src/io/parser/parser.hh
@@ -1,543 +1,544 @@
 /**
  * @file   parser.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Wed Nov 13 2013
  * @date last modification: Fri Apr 02 2021
  *
  * @brief  File parser interface
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2014-2021 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "aka_random_generator.hh"
 /* -------------------------------------------------------------------------- */
 #include <map>
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_PARSER_HH_
 #define AKANTU_PARSER_HH_
 
 namespace akantu {
 
 #if !defined(DOXYGEN)
 // clang-format off
 #define AKANTU_SECTION_TYPES                                            \
   (cohesive_inserter)                                                   \
   (contact)                                                             \
   (embedded_interface)                                                  \
   (friction)                                                            \
   (global)                                                              \
   (heat)                                                                \
   (integration_scheme)                                                  \
   (material)                                                            \
-  (phasefield)								\
+  (phasefield)                                                          \
+  (constitutive_law)                                                    \
   (mesh)                                                                \
   (model)                                                               \
   (model_solver)                                                        \
   (neighborhood)                                                        \
   (neighborhoods)                                                       \
   (non_linear_solver)                                                   \
   (non_local)                                                           \
   (rules)                                                               \
   (solver)                                                              \
   (time_step_solver)                                                    \
   (user)                                                                \
   (weight_function)                                                     \
   (contact_detector)							\
   (contact_resolution)							\
   (not_defined)
 // clang-format on
 
 /// Defines the possible section types
 AKANTU_CLASS_ENUM_DECLARE(ParserType, AKANTU_SECTION_TYPES)
 AKANTU_CLASS_ENUM_OUTPUT_STREAM(ParserType, AKANTU_SECTION_TYPES)
 AKANTU_CLASS_ENUM_INPUT_STREAM(ParserType, AKANTU_SECTION_TYPES)
 #else
 enum class ParserType {
   cohesive_inserter,
   contact,
   embedded_interface,
   friction,
   global,
   heat,
   integration_scheme,
   material,
   phasefield,
   mesh,
   model,
   model_solver,
   neighborhood,
   neighborhoods,
   non_linear_solver,
   non_local,
   rules,
   solver,
   time_step_solver,
   user,
   weight_function,
   not_defined
 };
 #endif
 
 /// Defines the possible search contexts/scopes (for parameter search)
 enum ParserParameterSearchCxt {
   _ppsc_current_scope = 0x1,
   _ppsc_parent_scope = 0x2,
   _ppsc_current_and_parent_scope = 0x3
 };
 /* ------------------------------------------------------------------------ */
 /* Parameters Class                                                         */
 /* ------------------------------------------------------------------------ */
 class ParserSection;
 
 /// @brief The ParserParameter objects represent the end of tree branches as
 /// they
 /// are the different informations contained in the input file.
 class ParserParameter {
 public:
   ParserParameter()
       : name(std::string()), value(std::string()), dbg_filename(std::string()) {
   }
 
   ParserParameter(const std::string & name, const std::string & value,
                   const ParserSection & parent_section)
       : parent_section(&parent_section), name(name), value(value),
         dbg_filename(std::string()) {}
 
   ParserParameter(const ParserParameter & param) = default;
 
   virtual ~ParserParameter() = default;
 
   /// Get parameter name
   const std::string & getName() const { return name; }
   /// Get parameter value
   const std::string & getValue() const { return value; }
 
   /// Set info for debug output
   void setDebugInfo(const std::string & filename, UInt line, UInt column) {
     dbg_filename = filename;
     dbg_line = line;
     dbg_column = column;
   }
 
   template <typename T> inline operator T() const;
 
   // template <typename T> inline operator Vector<T>() const;
   // template <typename T> inline operator Matrix<T>() const;
 
   /// Print parameter info in stream
   void printself(std::ostream & stream,
                  __attribute__((unused)) unsigned int indent = 0) const {
     stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line
            << ":" << dbg_column << ")";
   }
 
 private:
   void setParent(const ParserSection & sect) { parent_section = &sect; }
 
   friend class ParserSection;
 
 private:
   /// Pointer to the parent section
   const ParserSection * parent_section{nullptr};
   /// Name of the parameter
   std::string name;
   /// Value of the parameter
   std::string value;
   /// File for debug output
   std::string dbg_filename;
   /// Position of parameter in parsed file
   UInt dbg_line, dbg_column;
 };
 
 /* ------------------------------------------------------------------------ */
 /* Sections Class                                                           */
 /* ------------------------------------------------------------------------ */
 /// ParserSection represents a branch of the parsing tree.
 class ParserSection {
 public:
   using SubSections = std::multimap<ParserType, ParserSection>;
   using Parameters = std::map<std::string, ParserParameter>;
 
 private:
   using const_section_iterator_ = SubSections::const_iterator;
 
 public:
   /* ------------------------------------------------------------------------ */
   /* SubSection iterator                                                      */
   /* ------------------------------------------------------------------------ */
   /// Iterator on sections
   class const_section_iterator {
   public:
     using iterator_category = std::forward_iterator_tag;
     using value_type = ParserSection;
     using pointer = ParserSection *;
     using reference = ParserSection &;
 
     const_section_iterator() = default;
     const_section_iterator(const const_section_iterator_ & it) : it(it) {}
     const_section_iterator(const const_section_iterator & other) = default;
     const_section_iterator &
     operator=(const const_section_iterator & other) = default;
 
     const ParserSection & operator*() const { return it->second; }
     const ParserSection * operator->() const { return &(it->second); }
 
     bool operator==(const const_section_iterator & other) const {
       return it == other.it;
     }
 
     bool operator!=(const const_section_iterator & other) const {
       return it != other.it;
     }
 
     const_section_iterator & operator++() {
       ++it;
       return *this;
     }
 
     const_section_iterator operator++(int) {
       const_section_iterator tmp = *this;
       operator++();
       return tmp;
     }
 
   private:
     const_section_iterator_ it;
   };
 
   /* ------------------------------------------------------------------------ */
   /* Parameters iterator                                                      */
   /* ------------------------------------------------------------------------ */
   /// Iterator on parameters
   class const_parameter_iterator {
   public:
     const_parameter_iterator(const const_parameter_iterator & other) = default;
     const_parameter_iterator(const Parameters::const_iterator & it) : it(it) {}
 
     const_parameter_iterator &
     operator=(const const_parameter_iterator & other) {
       if (this != &other) {
         it = other.it;
       }
       return *this;
     }
 
     const ParserParameter & operator*() const { return it->second; }
     const ParserParameter * operator->() { return &(it->second); };
 
     bool operator==(const const_parameter_iterator & other) const {
       return it == other.it;
     }
 
     bool operator!=(const const_parameter_iterator & other) const {
       return it != other.it;
     }
 
     const_parameter_iterator & operator++() {
       ++it;
       return *this;
     }
 
     const_parameter_iterator operator++(int) {
       const_parameter_iterator tmp = *this;
       operator++();
       return tmp;
     }
 
   private:
     Parameters::const_iterator it;
   };
 
   /* ---------------------------------------------------------------------- */
   ParserSection() : name(std::string()) {}
 
   ParserSection(const std::string & name, ParserType type)
       : name(name), type(type) {}
 
   ParserSection(const std::string & name, ParserType type,
                 const std::string & option,
                 const ParserSection & parent_section)
       : parent_section(&parent_section), name(name), type(type),
         option(option) {}
 
   ParserSection(const ParserSection & section)
       : parent_section(section.parent_section), name(section.name),
         type(section.type), option(section.option),
         parameters(section.parameters),
         sub_sections_by_type(section.sub_sections_by_type) {
     setChldrenPointers();
   }
 
   ParserSection & operator=(const ParserSection & other) {
     if (&other != this) {
       parent_section = other.parent_section;
       name = other.name;
       type = other.type;
       option = other.option;
       parameters = other.parameters;
       sub_sections_by_type = other.sub_sections_by_type;
       setChldrenPointers();
     }
     return *this;
   }
 
   virtual ~ParserSection();
 
   virtual void printself(std::ostream & stream, unsigned int indent = 0) const;
 
   /* ---------------------------------------------------------------------- */
   /* Creation functions                                                     */
   /* ---------------------------------------------------------------------- */
 public:
   ParserParameter & addParameter(const ParserParameter & param);
   ParserSection & addSubSection(const ParserSection & section);
 
 protected:
   /// Clean ParserSection content
   void clean() {
     parameters.clear();
     sub_sections_by_type.clear();
   }
 
 private:
   void setChldrenPointers() {
     for (auto && param_pair : this->parameters) {
       param_pair.second.setParent(*this);
     }
 
     for (auto && sub_sect_pair : this->sub_sections_by_type) {
       sub_sect_pair.second.setParent(*this);
     }
   }
 
   /* ---------------------------------------------------------------------- */
   /* Accessors                                                              */
   /* ---------------------------------------------------------------------- */
 public:
   class SubSectionsRange
       : public std::pair<const_section_iterator, const_section_iterator> {
   public:
     SubSectionsRange(const const_section_iterator & first,
                      const const_section_iterator & second)
         : std::pair<const_section_iterator, const_section_iterator>(first,
                                                                     second) {}
     auto begin() { return this->first; }
     auto end() { return this->second; }
   };
 
   /// Get begin and end iterators on subsections of certain type
   auto getSubSections(ParserType type = ParserType::_not_defined) const {
     if (type != ParserType::_not_defined) {
       auto range = sub_sections_by_type.equal_range(type);
       return SubSectionsRange(range.first, range.second);
     }
     return SubSectionsRange(sub_sections_by_type.begin(),
                             sub_sections_by_type.end());
   }
 
   /// Get number of subsections of certain type
   UInt getNbSubSections(ParserType type = ParserType::_not_defined) const {
     if (type != ParserType::_not_defined) {
       return this->sub_sections_by_type.count(type);
     }
     return this->sub_sections_by_type.size();
   }
 
   /// Get begin and end iterators on parameters
   auto getParameters() const {
     return std::pair<const_parameter_iterator, const_parameter_iterator>(
         parameters.begin(), parameters.end());
   }
   /* ---------------------------------------------------------------------- */
   /// Get parameter within specified context
   const ParserParameter & getParameter(
       const std::string & name,
       ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const {
     Parameters::const_iterator it;
     if ((search_ctx & _ppsc_current_scope) != 0) {
       it = parameters.find(name);
     }
 
     if (it == parameters.end()) {
       if ((search_ctx & _ppsc_parent_scope) != 0 and
           parent_section != nullptr) {
         return parent_section->getParameter(name, search_ctx);
       }
 
       AKANTU_SILENT_EXCEPTION(
           "The parameter " << name
                            << " has not been found in the specified context");
     }
     return it->second;
   }
 
   /* ------------------------------------------------------------------------ */
   /// Get parameter within specified context, with a default value in case the
   /// parameter does not exists
   template <class T>
   T getParameter(
       const std::string & name, const T & default_value,
       ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const {
     try {
       T tmp = this->getParameter(name, search_ctx);
       return tmp;
     } catch (debug::Exception &) {
       return default_value;
     }
   }
 
   /* ------------------------------------------------------------------------ */
   /// Check if parameter exists within specified context
   bool hasParameter(
       const std::string & name,
       ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const {
 
     Parameters::const_iterator it;
     if ((search_ctx & _ppsc_current_scope) != 0) {
       it = parameters.find(name);
     }
 
     if (it == parameters.end()) {
       if ((search_ctx & _ppsc_parent_scope) != 0 and
           parent_section != nullptr) {
         return parent_section->hasParameter(name, search_ctx);
       }
       return false;
     }
     return true;
   }
 
   /* --------------------------------------------------------------------------
    */
   /// Get value of given parameter in context
   template <class T>
   T getParameterValue(
       const std::string & name,
       ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const {
     const ParserParameter & tmp_param = getParameter(name, search_ctx);
     T t = tmp_param;
     return t;
   }
 
   /* --------------------------------------------------------------------------
    */
   /// Get section name
   std::string getName() const { return name; }
   /// Get section type
   ParserType getType() const { return type; }
   /// Get section option
   std::string getOption(const std::string & def = "") const {
     return (not option.empty()) ? option : def;
   }
 
 protected:
   void setParent(const ParserSection & sect) { parent_section = &sect; }
 
   /* ---------------------------------------------------------------------- */
   /* Members                                                                */
   /* ---------------------------------------------------------------------- */
 private:
   /// Pointer to the parent section
   const ParserSection * parent_section{nullptr};
   /// Name of section
   std::string name;
   /// Type of section, see AKANTU_SECTION_TYPES
   ParserType type{ParserType::_not_defined};
   /// Section option
   std::string option;
   /// Map of parameters in section
   Parameters parameters;
   /// Multi-map of subsections
   SubSections sub_sections_by_type;
 };
 
 /* ------------------------------------------------------------------------ */
 /* Parser Class                                                             */
 /* ------------------------------------------------------------------------ */
 /// Root of parsing tree, represents the global ParserSection
 class Parser : public ParserSection {
 public:
   Parser() : ParserSection("global", ParserType::_global) {}
 
   void parse(const std::string & filename);
 
   std::string getLastParsedFile() const;
 
   static bool isPermissive() { return permissive_parser; }
 
 public:
   /// Parse real scalar
   static Real parseReal(const std::string & value,
                         const ParserSection & section);
   /// Parse real vector
   static Vector<Real> parseVector(const std::string & value,
                                   const ParserSection & section);
   /// Parse real matrix
   static Matrix<Real> parseMatrix(const std::string & value,
                                   const ParserSection & section);
   /// Parse real random parameter
   static RandomParameter<Real>
   parseRandomParameter(const std::string & value,
                        const ParserSection & section);
 
 protected:
   /// General parse function
   template <class T, class Grammar>
   static T parseType(const std::string & value, Grammar & grammar);
 
 protected:
   //  friend class Parsable;
   static bool permissive_parser;
   std::string last_parsed_file;
 };
 
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ParserParameter & _this) {
   _this.printself(stream);
   return stream;
 }
 
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ParserSection & section) {
   section.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 namespace std {
 template <> struct iterator_traits<::akantu::Parser::const_section_iterator> {
   using iterator_category = input_iterator_tag;
   using value_type = ::akantu::ParserParameter;
   using difference_type = ptrdiff_t;
   using pointer = const ::akantu::ParserParameter *;
   using reference = const ::akantu::ParserParameter &;
 };
 } // namespace std
 
 #include "parser_tmpl.hh"
 
 #endif /* AKANTU_PARSER_HH_ */
diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/common/constitutive_laws/constitutive_law_non_local_interface.hh
similarity index 67%
copy from src/model/solid_mechanics/materials/material_non_local.hh
copy to src/model/common/constitutive_laws/constitutive_law_non_local_interface.hh
index a5239c9b4..c58d22320 100644
--- a/src/model/solid_mechanics/materials/material_non_local.hh
+++ b/src/model/common/constitutive_laws/constitutive_law_non_local_interface.hh
@@ -1,119 +1,119 @@
 /**
- * @file   material_non_local.hh
+ * @file   constitutive_law_non_local.hh
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
- * @date last modification: Fri Apr 09 2021
+ * @date last modification: Mon Sep 11 2017
  *
- * @brief  Material class that handle the non locality of a law for example
- * damage.
+ * @brief  ConstitutiveLaw class that handle the non locality of a law for
+ * example damage.
  *
- *
- * @section LICENSE
- *
- * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Copyright (©)  2010-2018 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
+ * 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
+ * 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
+ * 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
+ * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
-#include "material.hh"
+#include "constitutive_law.hh"
 /* -------------------------------------------------------------------------- */
 
-#ifndef AKANTU_MATERIAL_NON_LOCAL_HH_
-#define AKANTU_MATERIAL_NON_LOCAL_HH_
+#ifndef AKANTU_CONSTITUTIVE_LAW_NON_LOCAL_HH_
+#define AKANTU_CONSTITUTIVE_LAW_NON_LOCAL_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
-class MaterialNonLocalInterface {
+class ConstitutiveLawNonLocalInterface {
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
-  /// initialize the material the non local parts of the material
-  void initMaterialNonLocal() {
+  /// initialize the constitutive_law the non local parts of the
+  /// constitutive_law
+  void initConstitutiveLawNonLocal() {
     this->registerNeighborhood();
     this->registerNonLocalVariables();
   };
 
   /// insert the quadrature points in the neighborhoods of the non-local manager
   virtual void insertIntegrationPointsInNeighborhoods(
       GhostType ghost_type,
       const ElementTypeMapReal & quadrature_points_coordinates) = 0;
 
   /// update the values in the non-local internal fields
   virtual void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened,
                                        const ID & field_id,
                                        GhostType ghost_type,
                                        ElementKind kind) = 0;
-  /// constitutive law
-  virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0;
 
 protected:
-  /// get the name of the neighborhood for this material
+  /// get the name of the neighborhood for this constitutive_law
   virtual ID getNeighborhoodName() = 0;
 
-  /// register the neighborhoods for the material
+  /// register the neighborhoods for the constitutive_law
   virtual void registerNeighborhood() = 0;
 
   /// register the non local internal variable
   virtual void registerNonLocalVariables() = 0;
 
   virtual inline void onElementsAdded(const Array<Element> & /*unused*/,
                                       const NewElementsEvent & /*unused*/) {}
 };
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
-template <UInt dim, class LocalParent>
-class MaterialNonLocal : public MaterialNonLocalInterface, public LocalParent {
+template <UInt dim, class ConstitutiveLawNonLocalInterface, class ConstitutiveLawParent>
+class ConstitutiveLawNonLocal : public ConstitutiveLawNonLocalInterface,
+                                public ConstitutiveLawParent {
+
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  explicit MaterialNonLocal(SolidMechanicsModel & model, const ID & id);
+  explicit ConstitutiveLawNonLocal(
+      typename ConstitutiveLawParent::ConstitutiveLawsHandler & handler,
+      const ID & id);
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// insert the quadrature points in the neighborhoods of the non-local manager
   void insertIntegrationPointsInNeighborhoods(
       GhostType ghost_type,
       const ElementTypeMapReal & quadrature_points_coordinates) override;
 
   /// update the values in the non-local internal fields
   void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened,
                                const ID & field_id, GhostType ghost_type,
                                ElementKind kind) override;
 
-  /// register the neighborhoods for the material
+  /// register the neighborhoods for the constitutive_law
   void registerNeighborhood() override;
 
 protected:
-  /// get the name of the neighborhood for this material
+  /// get the name of the neighborhood for this constitutive_law
   ID getNeighborhoodName() override { return this->name; }
 };
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
-#include "material_non_local_tmpl.hh"
+#include "constitutive_law_non_local_interface_tmpl.hh"
 
-#endif /* AKANTU_MATERIAL_NON_LOCAL_HH_ */
+#endif /* AKANTU_CONSTITUTIVE_LAW_NON_LOCAL_HH_ */
diff --git a/src/model/solid_mechanics/materials/material_non_local_tmpl.hh b/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
similarity index 77%
rename from src/model/solid_mechanics/materials/material_non_local_tmpl.hh
rename to src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
index 6c033166f..a28b6fa41 100644
--- a/src/model/solid_mechanics/materials/material_non_local_tmpl.hh
+++ b/src/model/common/constitutive_laws/constitutive_law_non_local_interface_tmpl.hh
@@ -1,129 +1,128 @@
 /**
- * @file   material_non_local_tmpl.hh
+ * @file   constitutive_law_non_local_tmpl.hh
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Jul 06 2017
  * @date last modification: Fri Mar 26 2021
  *
- * @brief  Implementation of material non-local
+ * @brief  Implementation of constitutive_law non-local
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2016-2021 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
-#include "material.hh"
-#include "material_non_local.hh"
+#include "constitutive_law.hh"
+#include "constitutive_law_non_local_interface.hh"
 #include "non_local_neighborhood.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
-template <UInt dim, class LocalParent>
-MaterialNonLocal<dim, LocalParent>::MaterialNonLocal(
-    SolidMechanicsModel & model, const ID & id)
-    : LocalParent(model, id) {
-  AKANTU_DEBUG_IN();
-
-  AKANTU_DEBUG_OUT();
+template <UInt dim, class ConstitutiveLawParent>
+ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::ConstitutiveLawNonLocal(
+    typename ConstitutiveLawParent::ConstitutiveLawsHandler & handler,
+    const ID & id)
+    : ConstitutiveLawParent(handler, id) {
 }
 
 /* -------------------------------------------------------------------------- */
-template <UInt dim, class LocalParent>
-void MaterialNonLocal<dim, LocalParent>::insertIntegrationPointsInNeighborhoods(
-    GhostType ghost_type,
-    const ElementTypeMapReal & quadrature_points_coordinates) {
+template <UInt dim, class ConstitutiveLawParent>
+void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::
+    insertIntegrationPointsInNeighborhoods(
+        GhostType ghost_type,
+        const ElementTypeMapReal & quadrature_points_coordinates) {
 
   IntegrationPoint q;
   q.ghost_type = ghost_type;
 
   auto & neighborhood = this->model.getNonLocalManager().getNeighborhood(
       this->getNeighborhoodName());
 
   for (auto & type :
        this->element_filter.elementTypes(dim, ghost_type, _ek_regular)) {
     q.type = type;
     const auto & elem_filter = this->element_filter(type, ghost_type);
     UInt nb_element = elem_filter.size();
 
     if (nb_element != 0U) {
       UInt nb_quad =
           this->getFEEngine().getNbIntegrationPoints(type, ghost_type);
 
       const auto & quads = quadrature_points_coordinates(type, ghost_type);
 
       auto nb_total_element =
           this->model.getMesh().getNbElement(type, ghost_type);
       auto quads_it = quads.begin_reinterpret(dim, nb_quad, nb_total_element);
       for (auto & elem : elem_filter) {
         Matrix<Real> quads = quads_it[elem];
         q.element = elem;
         for (UInt nq = 0; nq < nb_quad; ++nq) {
           q.num_point = nq;
           q.global_num = q.element * nb_quad + nq;
           neighborhood.insertIntegrationPoint(q, quads(nq));
         }
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
-template <UInt dim, class LocalParent>
-void MaterialNonLocal<dim, LocalParent>::updateNonLocalInternals(
+template <UInt dim, class ConstitutiveLawParent>
+void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::updateNonLocalInternals(
     ElementTypeMapReal & non_local_flattened, const ID & field_id,
     GhostType ghost_type, ElementKind kind) {
 
-  /// loop over all types in the material
+  /// loop over all types in the constitutive_law
   for (auto & el_type :
        this->element_filter.elementTypes(dim, ghost_type, kind)) {
     Array<Real> & internal =
         this->template getInternal<Real>(field_id)(el_type, ghost_type);
 
     auto & internal_flat = non_local_flattened(el_type, ghost_type);
     auto nb_component = internal_flat.getNbComponent();
 
     auto internal_it = internal.begin(nb_component);
     auto internal_flat_it = internal_flat.begin(nb_component);
 
     /// loop all elements for the given type
     const auto & filter = this->element_filter(el_type, ghost_type);
     UInt nb_quads =
         this->getFEEngine().getNbIntegrationPoints(el_type, ghost_type);
     for (auto & elem : filter) {
       for (UInt q = 0; q < nb_quads; ++q, ++internal_it) {
         UInt global_quad = elem * nb_quads + q;
         *internal_it = internal_flat_it[global_quad];
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
-template <UInt dim, class LocalParent>
-void MaterialNonLocal<dim, LocalParent>::registerNeighborhood() {
+template <UInt dim, class ConstitutiveLawParent>
+void ConstitutiveLawNonLocal<dim, ConstitutiveLawParent>::registerNeighborhood() {
   ID name = this->getNeighborhoodName();
-  this->model.getNonLocalManager().registerNeighborhood(name, name);
+  this->handler.getNonLocalManager().registerNeighborhood(name, name);
 }
 
 } // namespace akantu
diff --git a/src/model/common/constitutive_laws/constitutive_laws_handler.hh b/src/model/common/constitutive_laws/constitutive_laws_handler.hh
new file mode 100644
index 000000000..9c6b8ef4f
--- /dev/null
+++ b/src/model/common/constitutive_laws/constitutive_laws_handler.hh
@@ -0,0 +1,291 @@
+/**
+ * @file   constitutive_laws_handler.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Mohit Pundir <mohit.pundir@epfl.ch>
+ *
+ * @date creation  ven mar 26 2021
+ *
+ * @brief A handler for multiple consistutive laws
+ *
+ * @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/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#include "constitutive_law.hh"
+#include "mesh.hh"
+#include "mesh_events.hh"
+#include "non_local_manager_callback.hh"
+/* -------------------------------------------------------------------------- */
+
+#ifndef AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH
+#define AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH
+
+namespace akantu {
+template <class ConstitutiveLawsHandler> class ConstitutiveLawSelector;
+class NonLocalManager;
+} // namespace akantu
+
+/* -------------------------------------------------------------------------- */
+namespace akantu {
+
+template <class ConsistutiveLawsType>
+class ConstitutiveLawsHandler : public NonLocalManagerCallback {
+  /* ------------------------------------------------------------------------ */
+  /* Constructors/Destructors                                                 */
+  /* ------------------------------------------------------------------------ */
+public:
+  ConstitutiveLawsHandler(const Mesh & mesh, UInt spatial_dimension,
+                          const ID & parent_id)
+      : constitutive_law_index("constitutive_law index", id),
+        constitutive_law_local_numbering("constitutive_law local numbering",
+                                         id) {
+    constitutive_law_selector =
+        std::make_shared<DefaultConstitutiveLawSelector>(
+            constitutive_law_index);
+
+    constitutive_law_index.initialize(mesh, _element_kind = _ek_not_defined,
+                                      _default_value = UInt(-1),
+                                      _with_nb_element = true);
+    constitutive_law_local_numbering.initialize(
+        mesh, _element_kind = _ek_not_defined, _with_nb_element = true);
+  }
+
+  virtual ~ConstitutiveLawsHandler();
+
+  class NewConstitutiveLawElementsEvent : public NewElementsEvent {
+  public:
+    AKANTU_GET_MACRO_NOT_CONST(ConstitutiveLawList, constitutive_law,
+                               Array<UInt> &);
+    AKANTU_GET_MACRO(ConstitutiveLawList, constitutive_law,
+                     const Array<UInt> &);
+
+  protected:
+    Array<UInt> constitutive_law;
+  };
+
+  /* ------------------------------------------------------------------------ */
+  /* ConstitutiveLaws                                                         */
+  /* ------------------------------------------------------------------------ */
+public:
+  /// register an empty constitutive_law of a given type
+  auto & registerNewConstitutiveLaw(const ID & cl_name, const ID & cl_type,
+                                    const ID & opt_param);
+
+  /// reassigns constitutive_laws depending on the constitutive_law selector
+  virtual void reassignConstitutiveLaw();
+
+  template <class Func> void for_each_constitutive_law(Func && func) {
+    for (auto && constitutive_law : constitutive_laws) {
+      std::forward<Func>(func)(
+          aka::as_type<ConstitutiveLawType>(*constitutive_laws));
+    }
+  }
+
+protected:
+  /// register a constitutive_law in the dynamic database
+  auto & registerNewConstitutiveLaw(const ParserSection & cl_section);
+
+  /// read the constitutive_law files to instantiate all the constitutive_laws
+  void instantiateConstitutiveLaws();
+
+  /// set the element_id_by_constitutive_law and add the elements to the good
+  /// constitutive_laws
+  virtual void assignConstitutiveLawToElements(
+      const ElementTypeMapArray<UInt> * filter = nullptr);
+
+protected:
+  void splitElementByConstitutiveLaw(
+      const Array<Element> & elements,
+      std::vector<Array<Element>> & elements_per_cl) const;
+
+  template <typename Operation>
+  void splitByConstitutiveLaw(const Array<Element> & elements,
+                              Operation && op) const;
+
+  /* ------------------------------------------------------------------------ */
+  /* Dumpable interface (kept for convenience) and dumper relative functions  */
+  /* ------------------------------------------------------------------------ */
+public:
+  //! decide wether a field is a constitutive_law internal or not
+  bool isInternal(const std::string & field_name, ElementKind element_kind);
+
+  //! give the amount of data per element
+  virtual ElementTypeMap<UInt>
+  getInternalDataPerElem(const std::string & field_name, ElementKind kind);
+
+  //! flatten a given constitutive_law internal field
+  ElementTypeMapArray<Real> &
+  flattenInternal(const std::string & field_name, ElementKind kind,
+                  GhostType ghost_type = _not_ghost);
+
+  //! flatten all the registered constitutive_law internals
+  void flattenAllRegisteredInternals(ElementKind kind);
+
+  /* ------------------------------------------------------------------------ */
+  /* Methods                                                                  */
+  /* ------------------------------------------------------------------------ */
+public:
+  /* ------------------------------------------------------------------------ */
+  /* Accessors                                                                */
+  /* ------------------------------------------------------------------------ */
+public:
+  /// get an iterable on the constitutive_laws
+  inline decltype(auto) getConstitutiveLaws();
+
+  /// get an iterable on the constitutive_laws
+  inline decltype(auto) getConstitutiveLaws() const;
+
+  /// get a particular constitutive_law (by numerical constitutive_law index)
+  inline auto & getConstitutiveLaw(UInt cl_index);
+
+  /// get a particular constitutive_law (by numerical constitutive_law index)
+  inline const auto & getConstitutiveLaw(UInt cl_index) const;
+
+  /// get a particular constitutive_law (by constitutive_law name)
+  inline auto & getConstitutiveLaw(const std::string & name);
+
+  /// get a particular constitutive_law (by constitutive_law name)
+  inline const auto & getConstitutiveLaw(const std::string & name) const;
+
+  /// get a particular constitutive_law id from is name
+  inline UInt getConstitutiveLawIndex(const std::string & name) const;
+
+  /// give the number of constitutive_laws
+  inline UInt getNbConstitutiveLaws() const { return constitutive_laws.size(); }
+
+  /// give the constitutive_law internal index from its id
+  Int getInternalIndexFromID(const ID & id) const;
+
+  AKANTU_GET_MACRO(ConstitutiveLawByElement, constitutive_law_index,
+                   const ElementTypeMapArray<UInt> &);
+  AKANTU_GET_MACRO(ConstitutiveLawLocalNumbering,
+                   constitutive_law_local_numbering,
+                   const ElementTypeMapArray<UInt> &);
+
+  /// vectors containing local constitutive_law element index for each global
+  /// element index
+  AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConstitutiveLawByElement,
+                                         constitutive_law_index, UInt);
+  AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConstitutiveLawLocalNumbering,
+                                         constitutive_law_local_numbering,
+                                         UInt);
+
+  AKANTU_GET_MACRO_NOT_CONST(ConstitutiveLawSelector,
+                             *constitutive_law_selector,
+                             ConstitutiveLawSelector &);
+
+  /// Access the non_local_manager interface
+  AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &);
+
+  void setConstitutiveLawSelector(
+      std::shared_ptr<ConstitutiveLawSelector> constitutive_law_selector) {
+    this->constitutive_law_selector = std::move(constitutive_law_selector);
+  }
+
+  /* ------------------------------------------------------------------------ */
+  /* Data Accessor inherited members                                          */
+  /* ------------------------------------------------------------------------ */
+public:
+  UInt getNbData(const Array<Element> & elements,
+                 const SynchronizationTag & tag) const override;
+
+  void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
+                const SynchronizationTag & tag) const override;
+
+  void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
+                  const SynchronizationTag & tag) override;
+
+  /* ------------------------------------------------------------------------ */
+  /* Mesh Event Handler inherited members                                     */
+  /* ------------------------------------------------------------------------ */
+protected:
+  void onElementsAdded(const Array<Element> & element_list,
+                       const NewElementsEvent & event) override;
+  void onElementsRemoved(const Array<Element> & element_list,
+                         const ElementTypeMapArray<UInt> & new_numbering,
+                         const RemovedElementsEvent & event) override;
+
+  /* ------------------------------------------------------------------------ */
+  /* NonLocalManager inherited members                                        */
+  /* ------------------------------------------------------------------------ */
+protected:
+  void initializeNonLocal() override;
+
+  void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override;
+
+  void computeNonLocalStresses(GhostType ghost_type) override;
+
+  void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override;
+
+  /// update the values of the non local internal
+  void updateLocalInternal(ElementTypeMapReal & internal_flat,
+                           GhostType ghost_type, ElementKind kind) override;
+
+  /// copy the results of the averaging in the constitutive_laws
+  void updateNonLocalInternal(ElementTypeMapReal & internal_flat,
+                              GhostType ghost_type, ElementKind kind) override;
+
+  /* ------------------------------------------------------------------------ */
+  /* Class Members                                                            */
+  /* ------------------------------------------------------------------------ */
+private:
+  /// id of the derived class
+  const ID & parent_id;
+
+  /// underlying mesh
+  const Mesh & _mesh;
+
+  /// spatial_dimension of the derived model
+  UInt _spatial_dimension{_all_dimensions};
+
+  /// mapping between constitutive_law name and constitutive_law internal id
+  std::map<std::string, UInt> constitutive_laws_names_to_id;
+
+  /// Arrays containing the constitutive_law index for each element
+  ElementTypeMapArray<UInt> constitutive_law_index;
+
+  /// Arrays containing the position in the element filter of the
+  /// constitutive_law (constitutive_law's local numbering)
+  ElementTypeMapArray<UInt> constitutive_law_local_numbering;
+
+  /// list of used constitutive_laws
+  std::vector<std::unique_ptr<ConstitutiveLaw>> constitutive_laws;
+
+  /// class defining of to choose a constitutive_law
+  std::shared_ptr<ConstitutiveLawSelector> constitutive_law_selector;
+
+  using flatten_internal_map =
+      std::map<std::pair<std::string, ElementKind>,
+               std::unique_ptr<ElementTypeMapArray<Real>>>;
+
+  /// map a registered internals to be flattened for dump purposes
+  flatten_internal_map registered_internals;
+
+  /// tells if the constitutive_law are instantiated
+  bool are_constitutive_laws_instantiated{false};
+
+  /// non local manager
+  std::unique_ptr<NonLocalManager> non_local_manager;
+
+  template <class ConstitutiveLawsHandler> friend class ConstitutiveLaw;
+};
+
+} // namespace akantu
+#endif /* AKANTU_CONSTITUTIVE_LAWS_HANDLER_HH */
diff --git a/src/model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh b/src/model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh
new file mode 100644
index 000000000..b0b6052e9
--- /dev/null
+++ b/src/model/common/constitutive_laws/constitutive_laws_handler_tmpl.hh
@@ -0,0 +1,549 @@
+/**
+ * @file   constitutive_laws_handler_tmpl.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Mohit Pundir <mohit.pundir@epfl.ch>
+ *
+ * @date creation  ven mar 26 2021
+ *
+ * @brief Implementation of the consistutive laws handler.
+ *
+ * @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/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#include "constitutive_law_non_local_interface.hh"
+#include "constitutive_laws_handler.hh"
+#include "non_local_manager.hh"
+#include "parsable.hh"
+/* -------------------------------------------------------------------------- */
+
+#ifndef AKANTU_CONSTITUTIVE_LAWS_HANDLER_TMPL_HH
+#define AKANTU_CONSTITUTIVE_LAWS_HANDLER_TMPL_HH
+
+namespace akantu {
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline decltype(auto)
+ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaws() {
+  return make_transform_adaptor(
+      constitutive_laws, [](auto && value) -> decltype(auto) {
+        return aka::as_type<ConstitutiveLawType>(*value);
+      });
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline decltype(auto)
+ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaws() const {
+  return make_transform_adaptor(
+      constitutive_laws, [](auto && value) -> decltype(auto) {
+        return aka::as_type<ConstitutiveLawType>(*value);
+      });
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline UInt
+ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLawIndex(
+    const std::string & name) const {
+  auto it = constitutive_laws_names_to_id.find(name);
+  if (it == constitutive_laws_names_to_id.end()) {
+    AKANTU_SILENT_EXCEPTION(
+        "The model " << parent_id << " has no constitutive_law named " << name);
+  }
+
+  return it->second;
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline auto & ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw(
+    UInt cl_index) {
+  AKANTU_DEBUG_ASSERT(cl_index < constitutive_laws.size(),
+                      "The model " << parent_id
+                                   << " has no constitutive law number "
+                                   << cl_index);
+  return aka::as_type<ConstitutiveLawType>(*constitutive_laws.at(cl_index));
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline const auto &
+ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw(
+    UInt cl_index) const {
+  AKANTU_DEBUG_ASSERT(cl_index < constitutive_laws.size(),
+                      "The model " << parent_id
+                                   << " has no constitutive law number "
+                                   << cl_index);
+  return aka::as_type<ConstitutiveLawType>(*constitutive_laws.at(cl_index));
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline auto & ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw(
+    const std::string & name) {
+  std::map<std::string, UInt>::const_iterator it =
+      constitutive_laws_names_to_id.find(name);
+  if (it == constitutive_laws_names_to_id.end()) {
+    AKANTU_SILENT_EXCEPTION(
+        "The model " << parent_id << " has no constitutive_law named " << name);
+  }
+
+  return aka::as_type<ConstitutiveLawType>(*constitutive_laws[it->second]);
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+inline const auto &
+ConstitutiveLawsHandler<ConstitutiveLawType>::getConstitutiveLaw(
+    const std::string & name) const {
+  auto it = constitutive_laws_names_to_id.find(name);
+  if (it == constitutive_laws_names_to_id.end()) {
+    AKANTU_SILENT_EXCEPTION(
+        "The model " << parent_id << " has no constitutive_law named " << name);
+  }
+  return aka::as_type<ConstitutiveLawType>(*constitutive_laws[it->second]);
+}
+
+/* -------------------------------------------------------------------------- */
+#define AKA_FWD_CL(...) ::std::forward<decltype(__VA_ARGS__)>(__VA_ARGS__)
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+template <typename Operation>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::splitByConstitutiveLaw(
+    const Array<Element> & elements, Operation && op) const {
+  std::vector<Array<Element>> elements_per_cl(constitutive_laws.size());
+  this->splitElementByConstitutiveLaw(elements, elements_per_cl);
+
+  for (auto && cl : zip(constitutive_laws, elements_per_cl)) {
+    AKA_FWD_CL(op)(AKA_FWD_CL(*std::get<0>(cl)), AKA_FWD_CL(std::get<1>(cl)));
+  }
+}
+
+#undef AKA_FWD_CL
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+auto & ConstitutiveLawsHandler<ConstitutiveLawType>::registerNewConstitutiveLaw(
+    const ParserSection & section) {
+  std::string cl_name;
+  std::string cl_type = section.getName();
+  std::string opt_param = section.getOption();
+
+  try {
+    std::string tmp = section.getParameter("name");
+    cl_name = tmp; /** this can seam weird, but there is an ambiguous operator
+                    * overload that i couldn't solve. @todo remove the
+                    * weirdness of this code
+                    */
+  } catch (debug::Exception &) {
+    AKANTU_ERROR("A constitutive_law of type \'"
+                 << cl_type
+                 << "\' in the input file has been defined without a name!");
+  }
+
+  auto && cl = this->registerNewConstitutiveLaw(cl_name, cl_type, opt_param);
+
+  cl.parseSection(section);
+
+  return cl;
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+auto & ConstitutiveLawsHandler<ConstitutiveLawType>::registerNewConstitutiveLaw(
+    const ID & cl_name, const ID & cl_type, const ID & opt_param) {
+  AKANTU_DEBUG_ASSERT(constitutive_laws_names_to_id.find(cl_name) ==
+                          constitutive_laws_names_to_id.end(),
+                      "A constitutive_law with this name '"
+                          << cl_name << "' has already been registered. "
+                          << "Please use unique names for constitutive_laws");
+
+  UInt cl_count = constitutive_laws.size();
+  constitutive_laws_names_to_id[cl_name] = cl_count;
+
+  std::stringstream sstr_cl;
+  sstr_cl << this->id << ":" << cl_count << ":" << cl_type;
+  ID cl_id = sstr_cl.str();
+
+  std::unique_ptr<ConstitutiveLaw> constitutive_law =
+      ConstitutiveLawType::getFactory().allocate(cl_type, _spatial_dimension,
+                                                 opt_param, *this, cl_id);
+
+  constitutive_laws.push_back(std::move(constitutive_law));
+
+  return aka::as_type<ConstitutiveLawType>(*(constitutive_laws.back()));
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<
+    ConstitutiveLawType>::instantiateConstitutiveLaws() {
+  ParserSection model_section;
+  bool is_empty;
+  std::tie(model_section, is_empty) = this->getParserSection();
+
+  if (not is_empty) {
+    auto model_constitutive_laws =
+        model_section.getSubSections(ParserType::_constitutive_law);
+    for (const auto & section : model_constitutive_laws) {
+      this->registerNewConstitutiveLaw(section);
+    }
+  }
+
+  auto sub_sections =
+      this->parser.getSubSections(ParserType::_constitutive_law);
+  for (const auto & section : sub_sections) {
+    this->registerNewConstitutiveLaw(section);
+  }
+
+#ifdef AKANTU_DAMAGE_NON_LOCAL
+  for (auto & constitutive_law : constitutive_laws) {
+    if (dynamic_cast<ConstitutiveLawNonLocalInterface *>(
+            constitutive_law.get()) == nullptr) {
+      continue;
+    }
+
+    this->non_local_manager = std::make_unique<NonLocalManager>(
+        *this, *this, parent_id + ":non_local_manager");
+    break;
+  }
+#endif
+
+  if (constitutive_laws.empty()) {
+    AKANTU_EXCEPTION("No constitutive_laws where instantiated for the model"
+                     << parent_id);
+  }
+  are_constitutive_laws_instantiated = true;
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::
+    assignConstitutiveLawToElements(const ElementTypeMapArray<UInt> * filter) {
+
+  for_each_element(
+      _mesh,
+      [&](auto && element) {
+        UInt cl_index = (*constitutive_law_selector)(element);
+        AKANTU_DEBUG_ASSERT(cl_index < constitutive_laws.size(),
+                            "The constitutive_law selector returned an index "
+                            "that does not exists");
+        constitutive_law_index(element) = cl_index;
+      },
+      _element_filter = filter, _ghost_type = _not_ghost);
+
+  if (non_local_manager) {
+    non_local_manager->synchronize(*this,
+                                   SynchronizationTag::_constitutive_law_id);
+  }
+
+  for_each_element(
+      _mesh,
+      [&](auto && element) {
+        auto cl_index = constitutive_law_index(element);
+        auto index = constitutive_laws[cl_index]->addElement(element);
+        constitutive_law_local_numbering(element) = index;
+      },
+      _element_filter = filter, _ghost_type = _not_ghost);
+
+  // synchronize the element constitutive_law arrays
+  this->synchronize(SynchronizationTag::_constitutive_law_id);
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::initConstitutiveLaws() {
+  AKANTU_DEBUG_ASSERT(not constitutive_laws.empty(),
+                      "No constitutive_law to initialize !");
+
+  this->assignConstitutiveLawToElements();
+
+  for (auto & constitutive_law : constitutive_laws) {
+    /// init internals properties
+    constitutive_law->initConstitutiveLaw();
+  }
+
+  this->synchronize(SynchronizationTag::_smm_init_cl);
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+Int ConstitutiveLawsHandler<ConstitutiveLawType>::getInternalIndexFromID(
+    const ID & id) const {
+  AKANTU_DEBUG_IN();
+
+  auto it = constitutive_laws.begin();
+  auto end = constitutive_laws.end();
+
+  for (; it != end; ++it) {
+    if ((*it)->getID() == id) {
+      AKANTU_DEBUG_OUT();
+      return (it - constitutive_laws.begin());
+    }
+  }
+
+  AKANTU_DEBUG_OUT();
+  return -1;
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::reassignConstitutiveLaw() {
+  AKANTU_DEBUG_IN();
+
+  std::vector<Array<Element>> element_to_add(constitutive_laws.size());
+  std::vector<Array<Element>> element_to_remove(constitutive_laws.size());
+
+  for_each_element(mesh, [&](auto && element) {
+    auto old_constitutive_law = constitutive_law_index(element);
+    auto new_constitutive_law = (*constitutive_law_selector)(element);
+    if (old_constitutive_law != new_constitutive_law) {
+      element_to_add[new_constitutive_law].push_back(element);
+      element_to_remove[old_constitutive_law].push_back(element);
+    }
+  });
+
+  for (auto && data : enumerate(constitutive_laws)) {
+    auto cl_index = std::get<0>(data);
+    auto & cl = *std::get<1>(data);
+
+    cl.removeElements(element_to_remove[cl_index]);
+    cl.addElements(element_to_add[cl_index]);
+  }
+
+  AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::
+    updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) {
+  const ID field_name = criterion.getName();
+  for (auto & constitutive_law : constitutive_laws) {
+    if (!constitutive_law->isInternal<Real>(field_name, _ek_regular)) {
+      continue;
+    }
+
+    for (auto ghost_type : ghost_types) {
+      constitutive_law->flattenInternal(field_name, criterion, ghost_type,
+                                        _ek_regular);
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::initializeNonLocal() {
+  this->non_local_manager->synchronize(
+      *this, SynchronizationTag::_constitutive_law_id);
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::
+    insertIntegrationPointsInNeighborhoods(GhostType ghost_type) {
+  for (auto & cl : constitutive_laws) {
+    ConstitutiveLawNonLocalInterface * cl_non_local;
+    if ((cl_non_local = dynamic_cast<ConstitutiveLawNonLocalInterface *>(
+             cl.get())) == nullptr) {
+      continue;
+    }
+
+    ElementTypeMapArray<Real> quadrature_points_coordinates(
+        "quadrature_points_coordinates_tmp_nl", this->id);
+    quadrature_points_coordinates.initialize(this->getFEEngine(),
+                                             _nb_component = spatial_dimension,
+                                             _ghost_type = ghost_type);
+
+    for (const auto & type : quadrature_points_coordinates.elementTypes(
+             Model::spatial_dimension, ghost_type)) {
+      this->getFEEngine().computeIntegrationPointsCoordinates(
+          quadrature_points_coordinates(type, ghost_type), type, ghost_type);
+    }
+
+    cl_non_local->initConstitutiveLawNonLocal();
+
+    cl_non_local->insertIntegrationPointsInNeighborhoods(
+        ghost_type, quadrature_points_coordinates);
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::updateLocalInternal(
+    ElementTypeMapReal & internal_flat, GhostType ghost_type,
+    ElementKind kind) {
+  const ID field_name = internal_flat.getName();
+  for (auto & constitutive_law : constitutive_laws) {
+    if (constitutive_law->isInternal<Real>(field_name, kind)) {
+      constitutive_law->flattenInternal(field_name, internal_flat, ghost_type,
+                                        kind);
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::updateNonLocalInternal(
+    ElementTypeMapReal & internal_flat, GhostType ghost_type,
+    ElementKind kind) {
+
+  const ID field_name = internal_flat.getName();
+
+  for (auto & mat : constitutive_laws) {
+    if (not aka::is_of_type<ConstitutiveLawNonLocalInterface>(*mat)) {
+      continue;
+    }
+
+    auto & mat_non_local =
+        dynamic_cast<ConstitutiveLawNonLocalInterface &>(*mat);
+    mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
+                                          kind);
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::
+    splitElementByConstitutiveLaw(
+        const Array<Element> & elements,
+        std::vector<Array<Element>> & elements_per_cl) const {
+  for (const auto & el : elements) {
+    Element cl_el = el;
+    cl_el.element = this->constitutive_law_local_numbering(el);
+    elements_per_cl[this->constitutive_law_index(el)].push_back(cl_el);
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+UInt ConstitutiveLawsHandler<ConstitutiveLawType>::getNbData(
+    const Array<Element> & elements, const SynchronizationTag & tag) const {
+    UInt size{0};
+
+    if(tag == SynchronizationTag::_constitutive_law_id) {
+        size += elements.size() * sizeof(UInt);
+        break;
+    }
+
+
+  if (tag != SynchronizationTag::_constitutive_law_id) {
+    splitByConstitutiveLaw(elements, [&](auto && cl, auto && elements) {
+      size += cl.getNbData(elements, tag);
+    });
+  }
+  return size;
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::packData(
+    CommunicationBuffer & buffer, const Array<Element> & elements,
+    const SynchronizationTag & tag) const {
+  if (tag == SynchronizationTag::_constitutive_law_id) {
+    DataAccessor<Element>::packElementalDataHelper(
+        constitutive_law_index, buffer, elements, false, getFEEngine());
+  }
+
+  if (tag != SynchronizationTag::_constitutive_law_id) {
+    splitByConstitutiveLaw(elements, [&](auto && cl, auto && elements) {
+      cl.packData(buffer, elements, tag);
+    });
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::unpackData(
+    CommunicationBuffer & buffer, const Array<Element> & elements,
+    const SynchronizationTag & tag) {
+  if (tag == SynchronizationTag::_constitutive_law_id) {
+    for (auto && element : elements) {
+      UInt recv_cl_index;
+      buffer >> recv_cl_index;
+      UInt & cl_index = constitutive_law_index(element);
+      if (cl_index != UInt(-1)) {
+        continue;
+      }
+
+      // add ghosts element to the correct constitutive_law
+      cl_index = recv_cl_index;
+      UInt index = constitutive_laws[cl_index]->addElement(element);
+      constitutive_law_local_numbering(element) = index;
+    }
+  }
+
+  if (tag != SynchronizationTag::_constitutive_law_id) {
+    splitByConstitutiveLaw(elements, [&](auto && cl, auto && elements) {
+      cl.unpackData(buffer, elements, tag);
+    });
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::onElementsAdded(
+    const Array<Element> & element_list, const NewElementsEvent & event) {
+  this->constitutive_law_index.initialize(mesh, _element_kind = _ek_not_defined,
+                                          _with_nb_element = true,
+                                          _default_value = UInt(-1));
+  this->constitutive_law_local_numbering.initialize(
+      mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
+      _default_value = UInt(-1));
+
+  ElementTypeMapArray<UInt> filter("new_element_filter", this->getID());
+
+  for (const auto & elem : element_list) {
+    if (mesh.getSpatialDimension(elem.type) != spatial_dimension) {
+      continue;
+    }
+
+    if (!filter.exists(elem.type, elem.ghost_type)) {
+      filter.alloc(0, 1, elem.type, elem.ghost_type);
+    }
+    filter(elem.type, elem.ghost_type).push_back(elem.element);
+  }
+
+  // this fails in parallel if the event is sent on facet between constructor
+  // and initFull \todo: to debug...
+  this->assignConstitutiveLawToElements(&filter);
+
+  for (auto & constitutive_law : constitutive_laws) {
+    constitutive_law->onElementsAdded(element_list, event);
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+template <class ConstitutiveLawType>
+void ConstitutiveLawsHandler<ConstitutiveLawType>::onElementsRemoved(
+    const Array<Element> & element_list,
+    const ElementTypeMapArray<UInt> & new_numbering,
+    const RemovedElementsEvent & event) {
+  for (auto & constitutive_law : constitutive_laws) {
+    constitutive_law->onElementsRemoved(element_list, new_numbering, event);
+  }
+}
+
+} // namespace akantu
+
+#endif /* AKANTU_CONSTITUTIVE_LAWS_HANDLER_TMPL_HH */
diff --git a/src/model/common/consitutive_laws/internal_field.hh b/src/model/common/constitutive_laws/internal_field.hh
similarity index 100%
rename from src/model/common/consitutive_laws/internal_field.hh
rename to src/model/common/constitutive_laws/internal_field.hh
diff --git a/src/model/common/consitutive_laws/internal_field_tmpl.hh b/src/model/common/constitutive_laws/internal_field_tmpl.hh
similarity index 100%
rename from src/model/common/consitutive_laws/internal_field_tmpl.hh
rename to src/model/common/constitutive_laws/internal_field_tmpl.hh
diff --git a/src/model/common/consitutive_laws/random_internal_field.hh b/src/model/common/constitutive_laws/random_internal_field.hh
similarity index 100%
rename from src/model/common/consitutive_laws/random_internal_field.hh
rename to src/model/common/constitutive_laws/random_internal_field.hh
diff --git a/src/model/common/consitutive_laws/random_internal_field_tmpl.hh b/src/model/common/constitutive_laws/random_internal_field_tmpl.hh
similarity index 100%
rename from src/model/common/consitutive_laws/random_internal_field_tmpl.hh
rename to src/model/common/constitutive_laws/random_internal_field_tmpl.hh
diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/solid_mechanics/materials/material_non_local.hh
index a5239c9b4..cca6e20ec 100644
--- a/src/model/solid_mechanics/materials/material_non_local.hh
+++ b/src/model/solid_mechanics/materials/material_non_local.hh
@@ -1,119 +1,59 @@
 /**
  * @file   material_non_local.hh
  *
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Fri Apr 09 2021
  *
  * @brief  Material class that handle the non locality of a law for example
  * damage.
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
-#include "material.hh"
+#include "constitutive_law_non_local.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_MATERIAL_NON_LOCAL_HH_
 #define AKANTU_MATERIAL_NON_LOCAL_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
-class MaterialNonLocalInterface {
+class MaterialNonLocalInterface : public ConstitutiveLawNonLocalInterface {
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
-  /// initialize the material the non local parts of the material
-  void initMaterialNonLocal() {
-    this->registerNeighborhood();
-    this->registerNonLocalVariables();
-  };
-
-  /// insert the quadrature points in the neighborhoods of the non-local manager
-  virtual void insertIntegrationPointsInNeighborhoods(
-      GhostType ghost_type,
-      const ElementTypeMapReal & quadrature_points_coordinates) = 0;
-
-  /// update the values in the non-local internal fields
-  virtual void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened,
-                                       const ID & field_id,
-                                       GhostType ghost_type,
-                                       ElementKind kind) = 0;
   /// constitutive law
   virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0;
-
-protected:
-  /// get the name of the neighborhood for this material
-  virtual ID getNeighborhoodName() = 0;
-
-  /// register the neighborhoods for the material
-  virtual void registerNeighborhood() = 0;
-
-  /// register the non local internal variable
-  virtual void registerNonLocalVariables() = 0;
-
-  virtual inline void onElementsAdded(const Array<Element> & /*unused*/,
-                                      const NewElementsEvent & /*unused*/) {}
 };
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim, class LocalParent>
-class MaterialNonLocal : public MaterialNonLocalInterface, public LocalParent {
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-public:
-  explicit MaterialNonLocal(SolidMechanicsModel & model, const ID & id);
-
-  /* ------------------------------------------------------------------------ */
-  /* Methods                                                                  */
-  /* ------------------------------------------------------------------------ */
-public:
-  /// insert the quadrature points in the neighborhoods of the non-local manager
-  void insertIntegrationPointsInNeighborhoods(
-      GhostType ghost_type,
-      const ElementTypeMapReal & quadrature_points_coordinates) override;
-
-  /// update the values in the non-local internal fields
-  void updateNonLocalInternals(ElementTypeMapReal & non_local_flattened,
-                               const ID & field_id, GhostType ghost_type,
-                               ElementKind kind) override;
-
-  /// register the neighborhoods for the material
-  void registerNeighborhood() override;
-
-protected:
-  /// get the name of the neighborhood for this material
-  ID getNeighborhoodName() override { return this->name; }
-};
-
-} // namespace akantu
-
-/* -------------------------------------------------------------------------- */
-#include "material_non_local_tmpl.hh"
+using MaterialNonLocal =
+    ConstitutiveLawNonLocal<dim, MaterialNonLocalInterface, LocalParent>;
 
 #endif /* AKANTU_MATERIAL_NON_LOCAL_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 7136890bb..d83fd5211 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1259 +1,1108 @@
 /**
  * @file   solid_mechanics_model.cc
  *
  * @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Mauro Corrado <mauro.corrado@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Clement Roux <clement.roux@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Fri Apr 09 2021
  *
  * @brief  Implementation of the SolidMechanicsModel class
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "solid_mechanics_model.hh"
 #include "integrator_gauss.hh"
 #include "shape_lagrange.hh"
 #include "solid_mechanics_model_tmpl.hh"
 
 #include "communicator.hh"
 #include "element_synchronizer.hh"
 #include "sparse_matrix.hh"
 #include "synchronizer_registry.hh"
 
 #include "dumpable_inline_impl.hh"
 /* -------------------------------------------------------------------------- */
 #include "dumper_iohelper_paraview.hh"
 /* -------------------------------------------------------------------------- */
 #include "material_non_local.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /**
  * A solid mechanics model need a mesh  and a dimension to be created. the model
  * by it  self can not  do a lot,  the good init  functions should be  called in
  * order to configure the model depending on what we want to do.
  *
  * @param  mesh mesh  representing  the model  we  want to  simulate
  * @param dim spatial  dimension of the problem, if dim =  0 (default value) the
  * dimension of the problem is assumed to be the on of the mesh
  * @param id an id to identify the model
  * @param model_type this is an internal parameter for inheritance purposes
  */
-
 SolidMechanicsModel::SolidMechanicsModel(
     Mesh & mesh, UInt dim, const ID & id,
     std::shared_ptr<DOFManager> dof_manager, const ModelType model_type)
-    : Model(mesh, model_type, dim, id), material_index("material index", id),
-      material_local_numbering("material local numbering", id) {
+    : Model(mesh, model_type, dim, id) {
   AKANTU_DEBUG_IN();
 
   this->initDOFManager(dof_manager);
 
   this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
                                                Model::spatial_dimension);
 
   this->mesh.registerDumper<DumperParaview>("solid_mechanics_model", id, true);
   this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
                          _ek_regular);
 
-  material_selector = std::make_shared<DefaultMaterialSelector>(material_index);
+  // material_selector =
+  // std::make_shared<DefaultMaterialSelector>(material_index);
 
   this->registerDataAccessor(*this);
 
   if (this->mesh.isDistributed()) {
     auto & synchronizer = this->mesh.getElementSynchronizer();
     this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id);
     this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass);
     this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress);
     this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 SolidMechanicsModel::~SolidMechanicsModel() = default;
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) {
   Model::setTimeStep(time_step, solver_id);
 
   this->mesh.getDumper().setTimeStep(time_step);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Initialization                                                             */
 /* -------------------------------------------------------------------------- */
 /**
  * This function groups  many of the initialization in on  function. For most of
  * basics  case the  function should  be  enough. The  functions initialize  the
  * model, the internal  vectors, set them to 0, and  depending on the parameters
  * it also initialize the explicit or implicit solver.
  *
  * @param options
  * \parblock
  * contains the different options to initialize the model
  * \li \c analysis_method specify the type of solver to use
  * \endparblock
  */
 void SolidMechanicsModel::initFullImpl(const ModelOptions & options) {
-  material_index.initialize(mesh, _element_kind = _ek_not_defined,
-                            _default_value = UInt(-1), _with_nb_element = true);
-  material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
-                                      _with_nb_element = true);
-
   Model::initFullImpl(options);
 
   // initialize the materials
   if (not this->parser.getLastParsedFile().empty()) {
     this->instantiateMaterials();
     this->initMaterials();
   }
 
   this->initBC(*this, *displacement, *displacement_increment, *external_force);
 }
 
 /* -------------------------------------------------------------------------- */
 TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const {
   return TimeStepSolverType::_dynamic_lumped;
 }
 
 /* -------------------------------------------------------------------------- */
 ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
   case TimeStepSolverType::_dynamic_lumped: {
     options.non_linear_solver_type = NonLinearSolverType::_lumped;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_central_difference;
     options.solution_type["displacement"] = IntegrationScheme::_acceleration;
     break;
   }
   case TimeStepSolverType::_static: {
     options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
     options.integration_scheme_type["displacement"] =
         IntegrationSchemeType::_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
   case TimeStepSolverType::_dynamic: {
     if (this->method == _explicit_consistent_mass) {
       options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
       options.integration_scheme_type["displacement"] =
           IntegrationSchemeType::_central_difference;
       options.solution_type["displacement"] = IntegrationScheme::_acceleration;
     } else {
       options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
       options.integration_scheme_type["displacement"] =
           IntegrationSchemeType::_trapezoidal_rule_2;
       options.solution_type["displacement"] = IntegrationScheme::_displacement;
     }
     break;
   }
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 std::tuple<ID, TimeStepSolverType>
 SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
   switch (method) {
   case _explicit_lumped_mass: {
     return std::make_tuple("explicit_lumped",
                            TimeStepSolverType::_dynamic_lumped);
   }
   case _explicit_consistent_mass: {
     return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
   }
   case _static: {
     return std::make_tuple("static", TimeStepSolverType::_static);
   }
   case _implicit_dynamic: {
     return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
   }
   default:
     return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type,
                                      NonLinearSolverType /*unused*/) {
   auto & dof_manager = this->getDOFManager();
 
   /* ------------------------------------------------------------------------ */
   // for alloc type of solvers
   this->allocNodalField(this->displacement, spatial_dimension, "displacement");
   this->allocNodalField(this->previous_displacement, spatial_dimension,
                         "previous_displacement");
   this->allocNodalField(this->displacement_increment, spatial_dimension,
                         "displacement_increment");
   this->allocNodalField(this->internal_force, spatial_dimension,
                         "internal_force");
   this->allocNodalField(this->external_force, spatial_dimension,
                         "external_force");
   this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs");
   this->allocNodalField(this->current_position, spatial_dimension,
                         "current_position");
 
   // initialize the current positions
   this->current_position->copy(this->mesh.getNodes());
 
   /* ------------------------------------------------------------------------ */
   if (!dof_manager.hasDOFs("displacement")) {
     dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
     dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
     dof_manager.registerDOFsIncrement("displacement",
                                       *this->displacement_increment);
     dof_manager.registerDOFsPrevious("displacement",
                                      *this->previous_displacement);
   }
 
   /* ------------------------------------------------------------------------ */
   // for dynamic
   if (time_step_solver_type == TimeStepSolverType::_dynamic ||
       time_step_solver_type == TimeStepSolverType::_dynamic_lumped) {
     this->allocNodalField(this->velocity, spatial_dimension, "velocity");
     this->allocNodalField(this->acceleration, spatial_dimension,
                           "acceleration");
 
     if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
       dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
       dof_manager.registerDOFsDerivative("displacement", 2,
                                          *this->acceleration);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Initialize the model,basically it  pre-compute the shapes, shapes derivatives
  * and jacobian
  */
 void SolidMechanicsModel::initModel() {
   /// \todo add  the current position  as a parameter to  initShapeFunctions for
   /// large deformation
   getFEEngine().initShapeFunctions(_not_ghost);
   getFEEngine().initShapeFunctions(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleResidual() {
   AKANTU_DEBUG_IN();
 
   /* ------------------------------------------------------------------------ */
   // computes the internal forces
   this->assembleInternalForces();
 
   /* ------------------------------------------------------------------------ */
   this->getDOFManager().assembleToResidual("displacement",
                                            *this->external_force, 1);
   this->getDOFManager().assembleToResidual("displacement",
                                            *this->internal_force, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleResidual(const ID & residual_part) {
   AKANTU_DEBUG_IN();
 
   if ("external" == residual_part) {
     this->getDOFManager().assembleToResidual("displacement",
                                              *this->external_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   if ("internal" == residual_part) {
     this->assembleInternalForces();
     this->getDOFManager().assembleToResidual("displacement",
                                              *this->internal_force, 1);
     AKANTU_DEBUG_OUT();
     return;
   }
 
   AKANTU_CUSTOM_EXCEPTION(
       debug::SolverCallbackResidualPartUnknown(residual_part));
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) const {
   // \TODO check the materials to know what is the correct answer
   if (matrix_id == "C") {
     return _mt_not_defined;
   }
 
   if (matrix_id == "K") {
     auto matrix_type = _unsymmetric;
 
-    for (auto & material : materials) {
-      matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id));
-    }
+    for_each_constitutive_law([&](auto && material) {
+      matrix_type = std::max(matrix_type, material.getMatrixType(matrix_id));
+    });
   }
-  return _symmetric;
+}
+return _symmetric;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) {
   if (matrix_id == "K") {
     this->assembleStiffnessMatrix();
   } else if (matrix_id == "M") {
     this->assembleMass();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) {
   if (matrix_id == "M") {
     this->assembleMassLumped();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::beforeSolveStep() {
   for (auto & material : materials) {
     material->beforeSolveStep();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::afterSolveStep(bool converged) {
   for (auto & material : materials) {
     material->afterSolveStep(converged);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::predictor() { ++displacement_release; }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::corrector() { ++displacement_release; }
 
 /* -------------------------------------------------------------------------- */
 /**
  * This function computes the internal forces as \f$F_{int} = \int_{\Omega} N
  * \sigma d\Omega@\f$
  */
 void SolidMechanicsModel::assembleInternalForces() {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the internal forces");
 
   this->internal_force->zero();
 
   // compute the stresses of local elements
   AKANTU_DEBUG_INFO("Compute local stresses");
-  for (auto & material : materials) {
-    material->computeAllStresses(_not_ghost);
-  }
+  for_each_constitutive_law(
+      [](auto && material) { material.computeAllStresses(_not_ghost); });
 
   /* ------------------------------------------------------------------------ */
   /* Computation of the non local part */
   if (this->non_local_manager) {
     this->non_local_manager->computeAllNonLocalStresses();
   }
 
   // communicate the stresses
   AKANTU_DEBUG_INFO("Send data for residual assembly");
   this->asynchronousSynchronize(SynchronizationTag::_smm_stress);
 
   // assemble the forces due to local stresses
   AKANTU_DEBUG_INFO("Assemble residual for local elements");
-  for (auto & material : materials) {
-    material->assembleInternalForces(_not_ghost);
-  }
+  for_each_constitutive_law(
+      [](auto && material) { material.assembleInternalForces(_not_ghost); });
 
   // finalize communications
   AKANTU_DEBUG_INFO("Wait distant stresses");
   this->waitEndSynchronize(SynchronizationTag::_smm_stress);
 
   // assemble the stresses due to ghost elements
   AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
-  for (auto & material : materials) {
-    material->assembleInternalForces(_ghost);
-  }
+  for_each_constitutive_law(
+      [](auto && material) { material.assembleInternalForces(_ghost); });
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleStiffnessMatrix(bool need_to_reassemble) {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
 
   if (not this->getDOFManager().hasMatrix("K")) {
     this->getDOFManager().getNewMatrix("K", this->getMatrixType("K"));
   }
 
   // Check if materials need to recompute the matrix
-  for (auto & material : materials) {
+  for_each_constitutive_law([&](auto && material) {
     need_to_reassemble |= material->hasMatrixChanged("K");
-  }
+  });
 
   if (need_to_reassemble) {
     this->getDOFManager().getMatrix("K").zero();
 
     // call compute stiffness matrix on each local elements
-    for (auto & material : materials) {
-      material->assembleStiffnessMatrix(_not_ghost);
-    }
+    for_each_constitutive_law(
+        [](auto && material) { material.assembleStiffnessMatrix(_not_ghost); });
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateCurrentPosition() {
   if (this->current_position_release == this->displacement_release) {
     return;
   }
 
   this->current_position->copy(this->mesh.getNodes());
 
   auto cpos_it = this->current_position->begin(Model::spatial_dimension);
   auto cpos_end = this->current_position->end(Model::spatial_dimension);
   auto disp_it = this->displacement->begin(Model::spatial_dimension);
 
   for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
     *cpos_it += *disp_it;
   }
 
   this->current_position_release = this->displacement_release;
 }
 
 /* -------------------------------------------------------------------------- */
 const Array<Real> & SolidMechanicsModel::getCurrentPosition() {
   this->updateCurrentPosition();
   return *this->current_position;
 }
 
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::updateDataForNonLocalCriterion(
-    ElementTypeMapReal & criterion) {
-  const ID field_name = criterion.getName();
-  for (auto & material : materials) {
-    if (!material->isInternal<Real>(field_name, _ek_regular)) {
-      continue;
-    }
-
-    for (auto ghost_type : ghost_types) {
-      material->flattenInternal(field_name, criterion, ghost_type, _ek_regular);
-    }
-  }
-}
-
 /* -------------------------------------------------------------------------- */
 /* Information                                                                */
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getStableTimeStep() {
   AKANTU_DEBUG_IN();
 
   Real min_dt = getStableTimeStep(_not_ghost);
 
   /// reduction min over all processors
   mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min);
 
   AKANTU_DEBUG_OUT();
   return min_dt;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getStableTimeStep(GhostType ghost_type) {
   AKANTU_DEBUG_IN();
 
   Real min_dt = std::numeric_limits<Real>::max();
 
   this->updateCurrentPosition();
 
   Element elem;
   elem.ghost_type = ghost_type;
 
   for (auto type :
        mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) {
     elem.type = type;
     UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
     UInt nb_element = mesh.getNbElement(type);
 
-    auto mat_indexes = material_index(type, ghost_type).begin();
-    auto mat_loc_num = material_local_numbering(type, ghost_type).begin();
+    auto mat_indexes =
+        this->getConstitutiveLawByElement(type, ghost_type).begin();
+    auto mat_loc_num =
+        this->getConstitutiveLawLocalNumbering(type, ghost_type).begin();
 
     Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension);
     FEEngine::extractNodalToElementField(mesh, *current_position, X, type,
                                          _not_ghost);
 
     auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element);
 
     for (UInt el = 0; el < nb_element;
          ++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
       elem.element = *mat_loc_num;
       Real el_h = getFEEngine().getElementInradius(*X_el, type);
       Real el_c = this->materials[*mat_indexes]->getCelerity(elem);
       Real el_dt = el_h / el_c;
 
       min_dt = std::min(min_dt, el_dt);
     }
   }
 
   AKANTU_DEBUG_OUT();
   return min_dt;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getKineticEnergy() {
   AKANTU_DEBUG_IN();
 
   Real ekin = 0.;
   UInt nb_nodes = mesh.getNbNodes();
 
   if (this->getDOFManager().hasLumpedMatrix("M")) {
     this->assembleLumpedMatrix("M");
 
     auto m_it = this->mass->begin(Model::spatial_dimension);
     auto m_end = this->mass->end(Model::spatial_dimension);
     auto v_it = this->velocity->begin(Model::spatial_dimension);
 
     for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) {
       const auto & v = *v_it;
       const auto & m = *m_it;
 
       Real mv2 = 0.;
       auto is_local_node = mesh.isLocalOrMasterNode(n);
       // bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
       auto count_node = is_local_node; // && is_not_pbc_slave_node;
       if (count_node) {
         for (UInt i = 0; i < Model::spatial_dimension; ++i) {
           if (m(i) > std::numeric_limits<Real>::epsilon()) {
             mv2 += v(i) * v(i) * m(i);
           }
         }
       }
 
       ekin += mv2;
     }
   } else if (this->getDOFManager().hasMatrix("M")) {
     this->assembleMatrix("M");
 
     Array<Real> Mv(nb_nodes, Model::spatial_dimension);
     this->getDOFManager().assembleMatMulVectToArray("displacement", "M",
                                                     *this->velocity, Mv);
 
     for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension),
                             make_view(*this->velocity, spatial_dimension))) {
       ekin += std::get<2>(data).dot(std::get<1>(data)) *
               static_cast<Real>(mesh.isLocalOrMasterNode(std::get<0>(data)));
     }
   } else {
     AKANTU_ERROR("No function called to assemble the mass matrix.");
   }
 
   mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
 
   AKANTU_DEBUG_OUT();
   return ekin * .5;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getKineticEnergy(ElementType type, UInt index) {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
 
   Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension);
   Array<UInt> filter_element(1, 1, index);
 
   getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
                                                Model::spatial_dimension, type,
                                                _not_ghost, filter_element);
 
   auto vit = vel_on_quad.begin(Model::spatial_dimension);
   auto vend = vel_on_quad.end(Model::spatial_dimension);
 
   Vector<Real> rho_v2(nb_quadrature_points);
 
   Real rho = materials[material_index(type)(index)]->getRho();
 
   for (UInt q = 0; vit != vend; ++vit, ++q) {
     rho_v2(q) = rho * vit->dot(*vit);
   }
 
   AKANTU_DEBUG_OUT();
 
   return .5 * getFEEngine().integrate(rho_v2, type, index);
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getExternalWork() {
   AKANTU_DEBUG_IN();
 
   auto ext_force_it = external_force->begin(Model::spatial_dimension);
   auto int_force_it = internal_force->begin(Model::spatial_dimension);
   auto boun_it = blocked_dofs->begin(Model::spatial_dimension);
 
   decltype(ext_force_it) incr_or_velo_it;
   if (this->method == _static) {
     incr_or_velo_it =
         this->displacement_increment->begin(Model::spatial_dimension);
   } else {
     incr_or_velo_it = this->velocity->begin(Model::spatial_dimension);
   }
 
   Real work = 0.;
 
   UInt nb_nodes = this->mesh.getNbNodes();
 
   for (UInt n = 0; n < nb_nodes;
        ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) {
     const auto & int_force = *int_force_it;
     const auto & ext_force = *ext_force_it;
     const auto & boun = *boun_it;
     const auto & incr_or_velo = *incr_or_velo_it;
 
     bool is_local_node = this->mesh.isLocalOrMasterNode(n);
     // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
     bool count_node = is_local_node; // && is_not_pbc_slave_node;
 
     if (count_node) {
       for (UInt i = 0; i < Model::spatial_dimension; ++i) {
         if (boun(i)) {
           work -= int_force(i) * incr_or_velo(i);
         } else {
           work += ext_force(i) * incr_or_velo(i);
         }
       }
     }
   }
 
   mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum);
 
   if (this->method != _static) {
     work *= this->getTimeStep();
   }
 
   AKANTU_DEBUG_OUT();
   return work;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy();
   }
 
   if (energy_id == "external work") {
     return getExternalWork();
   }
 
   Real energy = 0.;
-  for (auto & material : materials) {
-    energy += material->getEnergy(energy_id);
-  }
+  for_each_constitutive_law(
+      [&](auto && material) { energy += material->getEnergy(energy_id); });
 
   /// reduction sum over all processors
   mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
 
   AKANTU_DEBUG_OUT();
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
                                     ElementType type, UInt index) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy(type, index);
   }
 
   UInt mat_index = this->material_index(type, _not_ghost)(index);
   UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
   Real energy =
       this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
 
   AKANTU_DEBUG_OUT();
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const ID & energy_id, const ID & group_id) {
   auto && group = mesh.getElementGroup(group_id);
   auto energy = 0.;
   for (auto && type : group.elementTypes()) {
     for (auto el : group.getElementsIterable(type)) {
       energy += getEnergy(energy_id, el);
     }
   }
 
   /// reduction sum over all processors
   mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
 
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
                                           const NewElementsEvent & event) {
-  AKANTU_DEBUG_IN();
-
-  this->material_index.initialize(mesh, _element_kind = _ek_not_defined,
-                                  _with_nb_element = true,
-                                  _default_value = UInt(-1));
-  this->material_local_numbering.initialize(
-      mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
-      _default_value = UInt(-1));
-
-  ElementTypeMapArray<UInt> filter("new_element_filter", this->getID());
-
-  for (const auto & elem : element_list) {
-    if (mesh.getSpatialDimension(elem.type) != spatial_dimension) {
-      continue;
-    }
-
-    if (!filter.exists(elem.type, elem.ghost_type)) {
-      filter.alloc(0, 1, elem.type, elem.ghost_type);
-    }
-    filter(elem.type, elem.ghost_type).push_back(elem.element);
-  }
-
-  // this fails in parallel if the event is sent on facet between constructor
-  // and initFull \todo: to debug...
-  this->assignMaterialToElements(&filter);
-
-  for (auto & material : materials) {
-    material->onElementsAdded(element_list, event);
-  }
-
-  AKANTU_DEBUG_OUT();
+  ConstitutiveLawsHandler<Material>::onElementsAdded(element_list, event);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onElementsRemoved(
     const Array<Element> & element_list,
     const ElementTypeMapArray<UInt> & new_numbering,
     const RemovedElementsEvent & event) {
-  for (auto & material : materials) {
-    material->onElementsRemoved(element_list, new_numbering, event);
-  }
+  ConstitutiveLawsHandler<Material>::onElementsRemoved(element_list,
+                                                       new_numbering, event);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
                                        const NewNodesEvent & event) {
   AKANTU_DEBUG_IN();
   UInt nb_nodes = mesh.getNbNodes();
 
   if (displacement) {
     displacement->resize(nb_nodes, 0.);
     ++displacement_release;
   }
   if (mass) {
     mass->resize(nb_nodes, 0.);
   }
   if (velocity) {
     velocity->resize(nb_nodes, 0.);
   }
   if (acceleration) {
     acceleration->resize(nb_nodes, 0.);
   }
   if (external_force) {
     external_force->resize(nb_nodes, 0.);
   }
   if (internal_force) {
     internal_force->resize(nb_nodes, 0.);
   }
   if (blocked_dofs) {
     blocked_dofs->resize(nb_nodes, false);
   }
   if (current_position) {
     current_position->resize(nb_nodes, 0.);
   }
 
   if (previous_displacement) {
     previous_displacement->resize(nb_nodes, 0.);
   }
   if (displacement_increment) {
     displacement_increment->resize(nb_nodes, 0.);
   }
 
   for (auto & material : materials) {
     material->onNodesAdded(nodes_list, event);
   }
 
   need_to_reassemble_lumped_mass = true;
   need_to_reassemble_mass = true;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/,
                                          const Array<UInt> & new_numbering,
                                          const RemovedNodesEvent & /*event*/) {
   if (displacement) {
     mesh.removeNodesFromArray(*displacement, new_numbering);
     ++displacement_release;
   }
   if (mass) {
     mesh.removeNodesFromArray(*mass, new_numbering);
   }
   if (velocity) {
     mesh.removeNodesFromArray(*velocity, new_numbering);
   }
   if (acceleration) {
     mesh.removeNodesFromArray(*acceleration, new_numbering);
   }
   if (internal_force) {
     mesh.removeNodesFromArray(*internal_force, new_numbering);
   }
   if (external_force) {
     mesh.removeNodesFromArray(*external_force, new_numbering);
   }
   if (blocked_dofs) {
     mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
   }
 
   // if (increment_acceleration)
   //   mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
   if (displacement_increment) {
     mesh.removeNodesFromArray(*displacement_increment, new_numbering);
   }
 
   if (previous_displacement) {
     mesh.removeNodesFromArray(*previous_displacement, new_numbering);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "Solid Mechanics Model [" << std::endl;
   stream << space << " + id                : " << id << std::endl;
   stream << space << " + spatial dimension : " << Model::spatial_dimension
          << std::endl;
 
   stream << space << " + fem [" << std::endl;
   getFEEngine().printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   stream << space << " + nodals information [" << std::endl;
   displacement->printself(stream, indent + 2);
   if (velocity) {
     velocity->printself(stream, indent + 2);
   }
   if (acceleration) {
     acceleration->printself(stream, indent + 2);
   }
   if (mass) {
     mass->printself(stream, indent + 2);
   }
   external_force->printself(stream, indent + 2);
   internal_force->printself(stream, indent + 2);
   blocked_dofs->printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   stream << space << " + material information [" << std::endl;
   material_index.printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   stream << space << " + materials [" << std::endl;
-  for (const auto & material : materials) {
-    material->printself(stream, indent + 2);
-  }
+  for_each_constitutive_law(
+      [&](auto && material) { material.printself(stream, indent + 2); });
   stream << space << " ]" << std::endl;
 
   stream << space << "]" << std::endl;
 }
 
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::initializeNonLocal() {
-  this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id);
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods(
-    GhostType ghost_type) {
-  for (auto & mat : materials) {
-    MaterialNonLocalInterface * mat_non_local;
-    if ((mat_non_local =
-             dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr) {
-      continue;
-    }
-
-    ElementTypeMapArray<Real> quadrature_points_coordinates(
-        "quadrature_points_coordinates_tmp_nl", this->id);
-    quadrature_points_coordinates.initialize(this->getFEEngine(),
-                                             _nb_component = spatial_dimension,
-                                             _ghost_type = ghost_type);
-
-    for (const auto & type : quadrature_points_coordinates.elementTypes(
-             Model::spatial_dimension, ghost_type)) {
-      this->getFEEngine().computeIntegrationPointsCoordinates(
-          quadrature_points_coordinates(type, ghost_type), type, ghost_type);
-    }
-
-    mat_non_local->initMaterialNonLocal();
-
-    mat_non_local->insertIntegrationPointsInNeighborhoods(
-        ghost_type, quadrature_points_coordinates);
-  }
-}
-
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::computeNonLocalStresses(GhostType ghost_type) {
-  for (auto & mat : materials) {
-    if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) {
+  for_each_constitutive_law([&](auto && material) {
+    if (not aka::is_of_type<MaterialNonLocalInterface>(mat)) {
       continue;
     }
 
-    auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
+    auto & mat_non_local = aka::as_type<MaterialNonLocalInterface>(mat);
     mat_non_local.computeNonLocalStresses(ghost_type);
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::updateLocalInternal(
-    ElementTypeMapReal & internal_flat, GhostType ghost_type,
-    ElementKind kind) {
-  const ID field_name = internal_flat.getName();
-  for (auto & material : materials) {
-    if (material->isInternal<Real>(field_name, kind)) {
-      material->flattenInternal(field_name, internal_flat, ghost_type, kind);
-    }
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::updateNonLocalInternal(
-    ElementTypeMapReal & internal_flat, GhostType ghost_type,
-    ElementKind kind) {
-
-  const ID field_name = internal_flat.getName();
-
-  for (auto & mat : materials) {
-    if (not aka::is_of_type<MaterialNonLocalInterface>(*mat)) {
-      continue;
-    }
-
-    auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
-    mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
-                                          kind);
-  }
+  });
 }
 
 /* -------------------------------------------------------------------------- */
 FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) {
   return getFEEngineClassBoundary<MyFEEngineType>(name);
 }
 
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::splitElementByMaterial(
-    const Array<Element> & elements,
-    std::vector<Array<Element>> & elements_per_mat) const {
-  for (const auto & el : elements) {
-    Element mat_el = el;
-    mat_el.element = this->material_local_numbering(el);
-    elements_per_mat[this->material_index(el)].push_back(mat_el);
-  }
-}
-
 /* -------------------------------------------------------------------------- */
 UInt SolidMechanicsModel::getNbData(const Array<Element> & elements,
                                     const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   UInt size = 0;
   UInt nb_nodes_per_element = 0;
 
   for (const Element & el : elements) {
     nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
   }
 
   switch (tag) {
-  case SynchronizationTag::_material_id: {
-    size += elements.size() * sizeof(UInt);
-    break;
-  }
   case SynchronizationTag::_smm_mass: {
     size += nb_nodes_per_element * sizeof(Real) *
             Model::spatial_dimension; // mass vector
     break;
   }
   case SynchronizationTag::_smm_for_gradu: {
     size += nb_nodes_per_element * Model::spatial_dimension *
             sizeof(Real); // displacement
     break;
   }
   case SynchronizationTag::_smm_boundary: {
     // force, displacement, boundary
     size += nb_nodes_per_element * Model::spatial_dimension *
             (2 * sizeof(Real) + sizeof(bool));
     break;
   }
   case SynchronizationTag::_for_dump: {
     // displacement, velocity, acceleration, residual, force
     size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5;
     break;
   }
   default: {
   }
   }
 
-  if (tag != SynchronizationTag::_material_id) {
-    splitByMaterial(elements, [&](auto && mat, auto && elements) {
-      size += mat.getNbData(elements, tag);
-    });
-  }
+  ConstitutiveLawsHandler<Material>::getNbData(elements, tag);
 
   AKANTU_DEBUG_OUT();
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
                                    const Array<Element> & elements,
                                    const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
-  case SynchronizationTag::_material_id: {
-    packElementalDataHelper(material_index, buffer, elements, false,
-                            getFEEngine());
-    break;
-  }
   case SynchronizationTag::_smm_mass: {
     packNodalDataHelper(*mass, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_for_gradu: {
     packNodalDataHelper(*displacement, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_for_dump: {
     packNodalDataHelper(*displacement, buffer, elements, mesh);
     packNodalDataHelper(*velocity, buffer, elements, mesh);
     packNodalDataHelper(*acceleration, buffer, elements, mesh);
     packNodalDataHelper(*internal_force, buffer, elements, mesh);
     packNodalDataHelper(*external_force, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_boundary: {
     packNodalDataHelper(*external_force, buffer, elements, mesh);
     packNodalDataHelper(*velocity, buffer, elements, mesh);
     packNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
     break;
   }
   default: {
   }
   }
 
-  if (tag != SynchronizationTag::_material_id) {
-    splitByMaterial(elements, [&](auto && mat, auto && elements) {
-      mat.packData(buffer, elements, tag);
-    });
-  }
+  ConstitutiveLawsHandler<Material>::packData(buffer, elements, tag);
+
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
                                      const Array<Element> & elements,
                                      const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
-  case SynchronizationTag::_material_id: {
-    for (auto && element : elements) {
-      UInt recv_mat_index;
-      buffer >> recv_mat_index;
-      UInt & mat_index = material_index(element);
-      if (mat_index != UInt(-1)) {
-        continue;
-      }
-
-      // add ghosts element to the correct material
-      mat_index = recv_mat_index;
-      UInt index = materials[mat_index]->addElement(element);
-      material_local_numbering(element) = index;
-    }
-    break;
-  }
   case SynchronizationTag::_smm_mass: {
     unpackNodalDataHelper(*mass, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_for_gradu: {
     unpackNodalDataHelper(*displacement, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_for_dump: {
     unpackNodalDataHelper(*displacement, buffer, elements, mesh);
     unpackNodalDataHelper(*velocity, buffer, elements, mesh);
     unpackNodalDataHelper(*acceleration, buffer, elements, mesh);
     unpackNodalDataHelper(*internal_force, buffer, elements, mesh);
     unpackNodalDataHelper(*external_force, buffer, elements, mesh);
     break;
   }
   case SynchronizationTag::_smm_boundary: {
     unpackNodalDataHelper(*external_force, buffer, elements, mesh);
     unpackNodalDataHelper(*velocity, buffer, elements, mesh);
     unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
     break;
   }
   default: {
   }
   }
 
-  if (tag != SynchronizationTag::_material_id) {
-    splitByMaterial(elements, [&](auto && mat, auto && elements) {
-      mat.unpackData(buffer, elements, tag);
-    });
-  }
-
+  ConstitutiveLawsHandler<Material>::unpackData(buffer, elements, tag);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs,
                                     const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   UInt size = 0;
   //  UInt nb_nodes = mesh.getNbNodes();
 
   switch (tag) {
   case SynchronizationTag::_smm_uv: {
     size += sizeof(Real) * Model::spatial_dimension * 2;
     break;
   }
   case SynchronizationTag::_smm_res: /* FALLTHRU */
   case SynchronizationTag::_smm_mass: {
     size += sizeof(Real) * Model::spatial_dimension;
     break;
   }
   case SynchronizationTag::_for_dump: {
     size += sizeof(Real) * Model::spatial_dimension * 5;
     break;
   }
   default: {
     AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
   return size * dofs.size();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
                                    const Array<UInt> & dofs,
                                    const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case SynchronizationTag::_smm_uv: {
     packDOFDataHelper(*displacement, buffer, dofs);
     packDOFDataHelper(*velocity, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_res: {
     packDOFDataHelper(*internal_force, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_mass: {
     packDOFDataHelper(*mass, buffer, dofs);
     break;
   }
   case SynchronizationTag::_for_dump: {
     packDOFDataHelper(*displacement, buffer, dofs);
     packDOFDataHelper(*velocity, buffer, dofs);
     packDOFDataHelper(*acceleration, buffer, dofs);
     packDOFDataHelper(*internal_force, buffer, dofs);
     packDOFDataHelper(*external_force, buffer, dofs);
     break;
   }
   default: {
     AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
                                      const Array<UInt> & dofs,
                                      const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case SynchronizationTag::_smm_uv: {
     unpackDOFDataHelper(*displacement, buffer, dofs);
     unpackDOFDataHelper(*velocity, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_res: {
     unpackDOFDataHelper(*internal_force, buffer, dofs);
     break;
   }
   case SynchronizationTag::_smm_mass: {
     unpackDOFDataHelper(*mass, buffer, dofs);
     break;
   }
   case SynchronizationTag::_for_dump: {
     unpackDOFDataHelper(*displacement, buffer, dofs);
     unpackDOFDataHelper(*velocity, buffer, dofs);
     unpackDOFDataHelper(*acceleration, buffer, dofs);
     unpackDOFDataHelper(*internal_force, buffer, dofs);
     unpackDOFDataHelper(*external_force, buffer, dofs);
     break;
   }
   default: {
     AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
+/* -------------------------------------------------------------------------- */
+void SolidMechanicsModel::applyEigenGradU(
+    const Matrix<Real> & prescribed_eigen_grad_u, const ID & material_name,
+    const GhostType ghost_type) {
+  AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() ==
+                          spatial_dimension * spatial_dimension,
+                      "The prescribed grad_u is not of the good size");
+  for (auto & material : this->getConstitutiveLaws()) {
+    if (material->getName() == material_name) {
+      material->applyEigenGradU(prescribed_eigen_grad_u, ghost_type);
+    }
+  }
+}
+
+/* -------------------------------------------------------------------------- */
+
 } // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index c59136979..290fee36b 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,601 +1,497 @@
 /**
  * @file   solid_mechanics_model.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Fri Apr 09 2021
  *
  * @brief  Model of Solid Mechanics
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "boundary_condition.hh"
+#include "constitutive_laws_handler.hh"
 #include "data_accessor.hh"
 #include "fe_engine.hh"
 #include "model.hh"
-#include "non_local_manager_callback.hh"
 #include "solid_mechanics_model_event_handler.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_
 #define AKANTU_SOLID_MECHANICS_MODEL_HH_
 
 namespace akantu {
 class Material;
 class MaterialSelector;
 class DumperIOHelper;
-class NonLocalManager;
 template <ElementKind kind, class IntegrationOrderFunctor>
 class IntegratorGauss;
 template <ElementKind kind> class ShapeLagrange;
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 class SolidMechanicsModel
     : public Model,
       public DataAccessor<Element>,
       public DataAccessor<UInt>,
       public BoundaryCondition<SolidMechanicsModel>,
-      public NonLocalManagerCallback,
+      public ConstitutiveLawsHandler<Material>,
       public EventHandlerManager<SolidMechanicsModelEventHandler> {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
-  class NewMaterialElementsEvent : public NewElementsEvent {
-  public:
-    AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &);
-    AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &);
-
-  protected:
-    Array<UInt> material;
-  };
-
   using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
 
 protected:
   using EventManager = EventHandlerManager<SolidMechanicsModelEventHandler>;
 
 public:
   SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions,
                       const ID & id = "solid_mechanics_model",
                       std::shared_ptr<DOFManager> dof_manager = nullptr,
                       ModelType model_type = ModelType::_solid_mechanics_model);
 
   ~SolidMechanicsModel() override;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 protected:
   /// initialize completely the model
   void initFullImpl(
       const ModelOptions & options = SolidMechanicsModelOptions()) override;
 
 public:
   /// initialize all internal arrays for materials
   virtual void initMaterials();
 
 protected:
   /// initialize the model
   void initModel() override;
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /// get some default values for derived classes
   std::tuple<ID, TimeStepSolverType>
   getDefaultSolverID(const AnalysisMethod & method) override;
 
   /* ------------------------------------------------------------------------ */
   /* Solver interface                                                         */
   /* ------------------------------------------------------------------------ */
 public:
   /// assembles the stiffness matrix,
   virtual void assembleStiffnessMatrix(bool need_to_reassemble = false);
   /// assembles the internal forces in the array internal_forces
   virtual void assembleInternalForces();
 
 protected:
   /// callback for the solver, this adds f_{ext} - f_{int} to the residual
   void assembleResidual() override;
 
   /// callback for the solver, this adds f_{ext} or  f_{int} to the residual
   void assembleResidual(const ID & residual_part) override;
   bool canSplitResidual() const override { return true; }
 
   /// get the type of matrix needed
   MatrixType getMatrixType(const ID & matrix_id) const override;
 
   /// callback for the solver, this assembles different matrices
   void assembleMatrix(const ID & matrix_id) override;
 
   /// callback for the solver, this assembles the stiffness matrix
   void assembleLumpedMatrix(const ID & matrix_id) override;
 
   /// callback for the solver, this is called at beginning of solve
   void predictor() override;
   /// callback for the solver, this is called at end of solve
   void corrector() override;
 
   /// callback for the solver, this is called at beginning of solve
   void beforeSolveStep() override;
   /// callback for the solver, this is called at end of solve
   void afterSolveStep(bool converged = true) override;
 
   /// Callback for the model to instantiate the matricees when needed
   void initSolver(TimeStepSolverType time_step_solver_type,
                   NonLinearSolverType non_linear_solver_type) override;
 
 protected:
   /* ------------------------------------------------------------------------ */
   TimeStepSolverType getDefaultSolverType() const override;
   /* ------------------------------------------------------------------------ */
   ModelSolverOptions
   getDefaultSolverOptions(const TimeStepSolverType & type) const override;
 
 public:
   bool isDefaultSolverExplicit() {
     return method == _explicit_lumped_mass ||
            method == _explicit_consistent_mass;
   }
 
 protected:
   /// update the current position vector
   void updateCurrentPosition();
 
   /* ------------------------------------------------------------------------ */
   /* Materials (solid_mechanics_model_material.cc)                            */
   /* ------------------------------------------------------------------------ */
 public:
-  /// register an empty material of a given type
-  Material & registerNewMaterial(const ID & mat_name, const ID & mat_type,
-                                 const ID & opt_param);
-
-  /// reassigns materials depending on the material selector
-  virtual void reassignMaterial();
-
   /// apply a constant eigen_grad_u on all quadrature points of a given material
   virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
                                const ID & material_name,
                                GhostType ghost_type = _not_ghost);
 
-protected:
-  /// register a material in the dynamic database
-  Material & registerNewMaterial(const ParserSection & mat_section);
-
-  /// read the material files to instantiate all the materials
-  void instantiateMaterials();
-
-  /// set the element_id_by_material and add the elements to the good materials
-  virtual void
-  assignMaterialToElements(const ElementTypeMapArray<UInt> * filter = nullptr);
-
   /* ------------------------------------------------------------------------ */
   /* Mass (solid_mechanics_model_mass.cc)                                     */
   /* ------------------------------------------------------------------------ */
 public:
   /// assemble the lumped mass matrix
   void assembleMassLumped();
 
   /// assemble the mass matrix for consistent mass resolutions
   void assembleMass();
 
 public:
   /// assemble the lumped mass matrix for local and ghost elements
   void assembleMassLumped(GhostType ghost_type);
 
   /// assemble the mass matrix for either _ghost or _not_ghost elements
   void assembleMass(GhostType ghost_type);
 
 protected:
   /// fill a vector of rho
   void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type);
 
   /// compute the kinetic energy
   Real getKineticEnergy();
   Real getKineticEnergy(ElementType type, UInt index);
 
   /// compute the external work (for impose displacement, the velocity should be
   /// given too)
   Real getExternalWork();
 
-  /* ------------------------------------------------------------------------ */
-  /* NonLocalManager inherited members                                        */
-  /* ------------------------------------------------------------------------ */
-protected:
-  void initializeNonLocal() override;
-
-  void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override;
-
-  void computeNonLocalStresses(GhostType ghost_type) override;
-
-  void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override;
-
-  /// update the values of the non local internal
-  void updateLocalInternal(ElementTypeMapReal & internal_flat,
-                           GhostType ghost_type, ElementKind kind) override;
-
-  /// copy the results of the averaging in the materials
-  void updateNonLocalInternal(ElementTypeMapReal & internal_flat,
-                              GhostType ghost_type, ElementKind kind) override;
-
   /* ------------------------------------------------------------------------ */
   /* Data Accessor inherited members                                          */
   /* ------------------------------------------------------------------------ */
 public:
   UInt getNbData(const Array<Element> & elements,
                  const SynchronizationTag & tag) const override;
 
   void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
                 const SynchronizationTag & tag) const override;
 
   void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
                   const SynchronizationTag & tag) override;
 
   UInt getNbData(const Array<UInt> & dofs,
                  const SynchronizationTag & tag) const override;
 
   void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
                 const SynchronizationTag & tag) const override;
 
   void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
                   const SynchronizationTag & tag) override;
 
-protected:
-  void
-  splitElementByMaterial(const Array<Element> & elements,
-                         std::vector<Array<Element>> & elements_per_mat) const;
-
-  template <typename Operation>
-  void splitByMaterial(const Array<Element> & elements, Operation && op) const;
-
   /* ------------------------------------------------------------------------ */
   /* Mesh Event Handler inherited members                                     */
   /* ------------------------------------------------------------------------ */
 protected:
   void onNodesAdded(const Array<UInt> & nodes_list,
                     const NewNodesEvent & event) override;
   void onNodesRemoved(const Array<UInt> & element_list,
                       const Array<UInt> & new_numbering,
                       const RemovedNodesEvent & event) override;
   void onElementsAdded(const Array<Element> & element_list,
                        const NewElementsEvent & event) override;
   void onElementsRemoved(const Array<Element> & element_list,
                          const ElementTypeMapArray<UInt> & new_numbering,
                          const RemovedElementsEvent & event) override;
   void onElementsChanged(const Array<Element> & /*unused*/,
                          const Array<Element> & /*unused*/,
                          const ElementTypeMapArray<UInt> & /*unused*/,
                          const ChangedElementsEvent & /*unused*/) override{};
 
   /* ------------------------------------------------------------------------ */
   /* Dumpable interface (kept for convenience) and dumper relative functions  */
   /* ------------------------------------------------------------------------ */
 public:
   virtual void onDump();
 
-  //! decide wether a field is a material internal or not
-  bool isInternal(const std::string & field_name, ElementKind element_kind);
-  //! give the amount of data per element
-  virtual ElementTypeMap<UInt>
-  getInternalDataPerElem(const std::string & field_name, ElementKind kind);
-
-  //! flatten a given material internal field
-  ElementTypeMapArray<Real> &
-  flattenInternal(const std::string & field_name, ElementKind kind,
-                  GhostType ghost_type = _not_ghost);
-  //! flatten all the registered material internals
-  void flattenAllRegisteredInternals(ElementKind kind);
-
-  //! inverse operation of the flatten
-  void inflateInternal(const std::string & field_name,
-                       const ElementTypeMapArray<Real> & field,
-                       ElementKind kind, GhostType ghost_type = _not_ghost);
-
   std::shared_ptr<dumpers::Field>
   createNodalFieldReal(const std::string & field_name,
                        const std::string & group_name,
                        bool padding_flag) override;
 
   std::shared_ptr<dumpers::Field>
   createNodalFieldBool(const std::string & field_name,
                        const std::string & group_name,
                        bool padding_flag) override;
 
   std::shared_ptr<dumpers::Field>
   createElementalField(const std::string & field_name,
                        const std::string & group_name, bool padding_flag,
                        UInt spatial_dimension, ElementKind kind) override;
 
   void dump(const std::string & dumper_name) override;
   void dump(const std::string & dumper_name, UInt step) override;
   void dump(const std::string & dumper_name, Real time, UInt step) override;
 
   void dump() override;
   void dump(UInt step) override;
   void dump(Real time, UInt step) override;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// set the value of the time step
   void setTimeStep(Real time_step, const ID & solver_id = "") override;
 
   /// get the value of the conversion from forces/ mass to acceleration
   AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
 
   /// set the value of the conversion from forces/ mass to acceleration
   AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
 
   /// get the SolidMechanicsModel::displacement array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement);
   /// get the SolidMechanicsModel::displacement array
   AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement);
 
   /// get the SolidMechanicsModel::previous_displacement array
   AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement);
 
   /// get the SolidMechanicsModel::current_position array
   const Array<Real> & getCurrentPosition();
 
   /// get  the SolidMechanicsModel::displacement_increment  array
   AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment);
   /// get  the SolidMechanicsModel::displacement_increment  array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment);
 
   /// get the lumped SolidMechanicsModel::mass array
   AKANTU_GET_MACRO_DEREF_PTR(Mass, mass);
 
   /// get the SolidMechanicsModel::velocity array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity);
   /// get the SolidMechanicsModel::velocity array
   AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity);
 
   /// get    the    SolidMechanicsModel::acceleration   array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration);
   /// get    the    SolidMechanicsModel::acceleration   array
   AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration);
 
   /// get the SolidMechanicsModel::external_force array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force);
   /// get the SolidMechanicsModel::external_force array
   AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force);
 
   /// get the SolidMechanicsModel::force array (external forces)
   [[deprecated("Use getExternalForce instead of this function")]] Array<Real> &
   getForce() {
     return getExternalForce();
   }
 
   /// get the SolidMechanicsModel::internal_force array (internal forces)
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force);
   /// get the SolidMechanicsModel::internal_force array (internal forces)
   AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force);
 
   /// get the SolidMechanicsModel::blocked_dofs array
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs);
   /// get the SolidMechanicsModel::blocked_dofs array
   AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs);
 
   /// get an iterable on the materials
   inline decltype(auto) getMaterials();
 
   /// get an iterable on the materials
   inline decltype(auto) getMaterials() const;
 
   /// get a particular material (by numerical material index)
   inline Material & getMaterial(UInt mat_index);
 
   /// get a particular material (by numerical material index)
   inline const Material & getMaterial(UInt mat_index) const;
 
   /// get a particular material (by material name)
   inline Material & getMaterial(const std::string & name);
 
   /// get a particular material (by material name)
   inline const Material & getMaterial(const std::string & name) const;
 
   /// get a particular material (by material name)
   inline const Material & getMaterial(const Element & element) const;
 
   /// get a particular material id from is name
   inline UInt getMaterialIndex(const std::string & name) const;
 
   /// give the number of materials
-  inline UInt getNbMaterials() const { return materials.size(); }
-
-  /// give the material internal index from its id
-  Int getInternalIndexFromID(const ID & id) const;
+  inline UInt getNbMaterials() const { return getNbConstitutiveLaws(); }
 
   /// compute the stable time step
   Real getStableTimeStep();
 
   /**
    * @brief Returns the total energy for a given energy type
    *
    * Energy types of SolidMechanicsModel expected as argument are:
    *   - `kinetic`
    *   - `external work`
    *
    * Other energy types are passed on to the materials. All materials should
    * define a `potential` energy type. For additional energy types, see material
    * documentation.
    */
   Real getEnergy(const std::string & energy_id);
 
   /// Compute energy for an element type and material index
   Real getEnergy(const std::string & energy_id, ElementType type, UInt index);
 
   /// Compute energy for an individual element
   Real getEnergy(const std::string & energy_id, const Element & element) {
     return getEnergy(energy_id, element.type, element.element);
   }
 
   /// Compute energy for an element group
   Real getEnergy(const ID & energy_id, const ID & group_id);
 
-  AKANTU_GET_MACRO(MaterialByElement, material_index,
-                   const ElementTypeMapArray<UInt> &);
-  AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering,
-                   const ElementTypeMapArray<UInt> &);
-
-  /// vectors containing local material element index for each global element
-  /// index
-  AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index,
-                                         UInt);
-  // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt);
-  AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering,
-                                         material_local_numbering, UInt);
-  // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering,
-  //                                  material_local_numbering, UInt);
-
-  AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, material_selector,
-                             std::shared_ptr<MaterialSelector>);
+  // this function is kept for backward compatinility
+  decltype(auto) getMaterialByElement() const {
+    return this->getConstitutiveLawByElement();
+  }
+
+  // this function is kept for backward compatinility
+  decltype(auto) getMaterialLocalNumbering() const {
+    return this->getConstitutiveLawLocalNumbering();
+  }
+
+  // this function is kept for backward compatinility
+  decltype(auto) getMaterialByElement(ElementType type,
+                                      GhostType ghost_type = _not_ghost) const {
+    return this->getConstitutiveLawByElement(type, ghost_type);
+  }
+
+  // this function is kept for backward compatinility
+  decltype(auto)
+  getMaterialLocalNumbering(ElementType type,
+                            GhostType ghost_type = _not_ghost) const {
+    return this->getConstitutiveLawLocalNumbering(type, ghost_type);
+  }
+
+  // this function is kept for backward compatinility
+  decltype(auto) getMaterialSelector() {
+    return this->getConstitutiveLawSelector();
+  }
+
+  // this function is kept for backward compatinility
   void
   setMaterialSelector(std::shared_ptr<MaterialSelector> material_selector) {
-    this->material_selector = std::move(material_selector);
+    this->setConstitutiveLawSelector(std::move(material_selector));
   }
 
-  /// Access the non_local_manager interface
-  AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &);
-
   /// get the FEEngine object to integrate or interpolate on the boundary
   FEEngine & getFEEngineBoundary(const ID & name = "") override;
 
 protected:
   /// compute the stable time step
   Real getStableTimeStep(GhostType ghost_type);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   /// release version of the displacement array
   UInt displacement_release{0};
 
   /// release version of the current_position array
   UInt current_position_release{0};
 
   /// Check if materials need to recompute the mass array
   bool need_to_reassemble_lumped_mass{true};
   /// Check if materials need to recompute the mass matrix
   bool need_to_reassemble_mass{true};
 
-  /// mapping between material name and material internal id
-  std::map<std::string, UInt> materials_names_to_id;
-
 protected:
   /// conversion coefficient form force/mass to acceleration
   Real f_m2a{1.0};
 
   /// displacements array
   std::unique_ptr<Array<Real>> displacement;
 
   /// displacements array at the previous time step (used in finite deformation)
   std::unique_ptr<Array<Real>> previous_displacement;
 
   /// increment of displacement
   std::unique_ptr<Array<Real>> displacement_increment;
 
   /// lumped mass array
   std::unique_ptr<Array<Real>> mass;
 
   /// velocities array
   std::unique_ptr<Array<Real>> velocity;
 
   /// accelerations array
   std::unique_ptr<Array<Real>> acceleration;
 
   /// external forces array
   std::unique_ptr<Array<Real>> external_force;
 
   /// internal forces array
   std::unique_ptr<Array<Real>> internal_force;
 
   /// array specifing if a degree of freedom is blocked or not
   std::unique_ptr<Array<bool>> blocked_dofs;
 
   /// array of current position used during update residual
   std::unique_ptr<Array<Real>> current_position;
-
-  /// Arrays containing the material index for each element
-  ElementTypeMapArray<UInt> material_index;
-
-  /// Arrays containing the position in the element filter of the material
-  /// (material's local numbering)
-  ElementTypeMapArray<UInt> material_local_numbering;
-
-  /// list of used materials
-  std::vector<std::unique_ptr<Material>> materials;
-
-  /// class defining of to choose a material
-  std::shared_ptr<MaterialSelector> material_selector;
-
-  using flatten_internal_map =
-      std::map<std::pair<std::string, ElementKind>,
-               std::unique_ptr<ElementTypeMapArray<Real>>>;
-
-  /// tells if the material are instantiated
-  flatten_internal_map registered_internals;
-
-  /// non local manager
-  std::unique_ptr<NonLocalManager> non_local_manager;
-
-  /// tells if the material are instantiated
-  bool are_materials_instantiated{false};
-
-  friend class Material;
-
-  template <class Model_> friend class CouplerSolidContactTemplate;
 };
 
 /* -------------------------------------------------------------------------- */
 namespace BC {
   namespace Neumann {
     using FromStress = FromHigherDim;
     using FromTraction = FromSameDim;
   } // namespace Neumann
 } // namespace BC
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "material.hh"
 #include "parser.hh"
-
-#include "solid_mechanics_model_inline_impl.hh"
-#include "solid_mechanics_model_tmpl.hh"
 /* -------------------------------------------------------------------------- */
 
 #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.hh
deleted file mode 100644
index d3dfed74f..000000000
--- a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.hh
+++ /dev/null
@@ -1,118 +0,0 @@
-/**
- * @file   solid_mechanics_model_inline_impl.hh
- *
- * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
- * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
- *
- * @date creation: Wed Aug 04 2010
- * @date last modification: Fri Mar 26 2021
- *
- * @brief  Implementation of the inline functions of the SolidMechanicsModel
- * class
- *
- *
- * @section LICENSE
- *
- * Copyright (©) 2015-2021 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/>.
- *
- */
-
-/* -------------------------------------------------------------------------- */
-#include "aka_named_argument.hh"
-#include "material_selector.hh"
-#include "material_selector_tmpl.hh"
-#include "solid_mechanics_model.hh"
-/* -------------------------------------------------------------------------- */
-
-#ifndef AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_HH_
-#define AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_HH_
-
-namespace akantu {
-
-/* -------------------------------------------------------------------------- */
-inline decltype(auto) SolidMechanicsModel::getMaterials() {
-  return make_dereference_adaptor(materials);
-}
-
-/* -------------------------------------------------------------------------- */
-inline decltype(auto) SolidMechanicsModel::getMaterials() const {
-  return make_dereference_adaptor(materials);
-}
-
-/* -------------------------------------------------------------------------- */
-inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) {
-  AKANTU_DEBUG_ASSERT(mat_index < materials.size(),
-                      "The model " << id << " has no material no "
-                                   << mat_index);
-  return *materials.at(mat_index);
-}
-
-/* -------------------------------------------------------------------------- */
-inline const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const {
-  AKANTU_DEBUG_ASSERT(mat_index < materials.size(),
-                      "The model " << id << " has no material no "
-                                   << mat_index);
-  return *materials.at(mat_index);
-}
-
-/* -------------------------------------------------------------------------- */
-inline Material & SolidMechanicsModel::getMaterial(const std::string & name) {
-  std::map<std::string, UInt>::const_iterator it =
-      materials_names_to_id.find(name);
-  if (it == materials_names_to_id.end()) {
-    AKANTU_SILENT_EXCEPTION("The model " << id << " has no material named "
-                                         << name);
-  }
-
-  return *materials[it->second];
-}
-
-/* -------------------------------------------------------------------------- */
-inline const Material &
-SolidMechanicsModel::getMaterial(const Element & element) const {
-  auto mat_id = material_index(element);
-  return *materials[mat_id];
-}
-
-/* -------------------------------------------------------------------------- */
-inline UInt
-SolidMechanicsModel::getMaterialIndex(const std::string & name) const {
-  auto it = materials_names_to_id.find(name);
-  if (it == materials_names_to_id.end()) {
-    AKANTU_SILENT_EXCEPTION("The model " << id << " has no material named "
-                                         << name);
-  }
-
-  return it->second;
-}
-
-/* -------------------------------------------------------------------------- */
-inline const Material &
-SolidMechanicsModel::getMaterial(const std::string & name) const {
-  auto it = materials_names_to_id.find(name);
-  if (it == materials_names_to_id.end()) {
-    AKANTU_SILENT_EXCEPTION("The model " << id << " has no material named "
-                                         << name);
-  }
-  return *materials[it->second];
-}
-
-/* -------------------------------------------------------------------------- */
-} // namespace akantu
-
-#endif /* AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_HH_ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_material.cc b/src/model/solid_mechanics/solid_mechanics_model_material.cc
deleted file mode 100644
index c06fea251..000000000
--- a/src/model/solid_mechanics/solid_mechanics_model_material.cc
+++ /dev/null
@@ -1,247 +0,0 @@
-/**
- * @file   solid_mechanics_model_material.cc
- *
- * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
- * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
- *
- * @date creation: Fri Nov 26 2010
- * @date last modification: Fri Mar 26 2021
- *
- * @brief  instatiation of materials
- *
- *
- * @section LICENSE
- *
- * Copyright (©) 2010-2021 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/>.
- *
- */
-
-/* -------------------------------------------------------------------------- */
-#include "aka_factory.hh"
-#include "aka_math.hh"
-#include "material_non_local.hh"
-#include "mesh_iterators.hh"
-#include "non_local_manager.hh"
-#include "solid_mechanics_model.hh"
-/* -------------------------------------------------------------------------- */
-
-namespace akantu {
-
-/* -------------------------------------------------------------------------- */
-Material &
-SolidMechanicsModel::registerNewMaterial(const ParserSection & section) {
-  std::string mat_name;
-  std::string mat_type = section.getName();
-  std::string opt_param = section.getOption();
-
-  try {
-    std::string tmp = section.getParameter("name");
-    mat_name = tmp; /** this can seam weird, but there is an ambiguous operator
-                     * overload that i couldn't solve. @todo remove the
-                     * weirdness of this code
-                     */
-  } catch (debug::Exception &) {
-    AKANTU_ERROR("A material of type \'"
-                 << mat_type
-                 << "\' in the input file has been defined without a name!");
-  }
-  Material & mat = this->registerNewMaterial(mat_name, mat_type, opt_param);
-
-  mat.parseSection(section);
-
-  return mat;
-}
-
-/* -------------------------------------------------------------------------- */
-Material & SolidMechanicsModel::registerNewMaterial(const ID & mat_name,
-                                                    const ID & mat_type,
-                                                    const ID & opt_param) {
-  AKANTU_DEBUG_ASSERT(materials_names_to_id.find(mat_name) ==
-                          materials_names_to_id.end(),
-                      "A material with this name '"
-                          << mat_name << "' has already been registered. "
-                          << "Please use unique names for materials");
-
-  UInt mat_count = materials.size();
-  materials_names_to_id[mat_name] = mat_count;
-
-  ID mat_id = this->id + ":" + std::to_string(mat_count) + ":" + mat_type;
-
-  std::unique_ptr<Material> material = MaterialFactory::getInstance().allocate(
-      mat_type, spatial_dimension, opt_param, *this, mat_id);
-  material->setName(mat_name);
-  materials.push_back(std::move(material));
-
-  return *(materials.back());
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::instantiateMaterials() {
-  ParserSection model_section;
-  bool is_empty;
-  std::tie(model_section, is_empty) = this->getParserSection();
-
-  if (not is_empty) {
-    auto model_materials = model_section.getSubSections(ParserType::_material);
-    for (const auto & section : model_materials) {
-      this->registerNewMaterial(section);
-    }
-  }
-
-  auto sub_sections = this->parser.getSubSections(ParserType::_material);
-  for (const auto & section : sub_sections) {
-    this->registerNewMaterial(section);
-  }
-
-#ifdef AKANTU_DAMAGE_NON_LOCAL
-  for (auto & material : materials) {
-    if (dynamic_cast<MaterialNonLocalInterface *>(material.get()) == nullptr) {
-      continue;
-    }
-
-    this->non_local_manager = std::make_unique<NonLocalManager>(
-        *this, *this, id + ":non_local_manager");
-    break;
-  }
-#endif
-
-  if (materials.empty()) {
-    AKANTU_EXCEPTION("No materials where instantiated for the model"
-                     << getID());
-  }
-  are_materials_instantiated = true;
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::assignMaterialToElements(
-    const ElementTypeMapArray<UInt> * filter) {
-
-  for_each_element(
-      mesh,
-      [&](auto && element) {
-        UInt mat_index = (*material_selector)(element);
-        AKANTU_DEBUG_ASSERT(
-            mat_index < materials.size(),
-            "The material selector returned an index that does not exists");
-        material_index(element) = mat_index;
-      },
-      _element_filter = filter, _ghost_type = _not_ghost);
-
-  if (non_local_manager) {
-    non_local_manager->synchronize(*this, SynchronizationTag::_material_id);
-  }
-
-  for_each_element(
-      mesh,
-      [&](auto && element) {
-        auto mat_index = material_index(element);
-        auto index = materials[mat_index]->addElement(element);
-        material_local_numbering(element) = index;
-      },
-      _element_filter = filter, _ghost_type = _not_ghost);
-
-  // synchronize the element material arrays
-  this->synchronize(SynchronizationTag::_material_id);
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::initMaterials() {
-  AKANTU_DEBUG_ASSERT(not materials.empty(), "No material to initialize !");
-
-  // if (!are_materials_instantiated)
-  //   instantiateMaterials();
-
-  this->assignMaterialToElements();
-
-  for (auto & material : materials) {
-    /// init internals properties
-    material->initMaterial();
-  }
-
-  this->synchronize(SynchronizationTag::_smm_init_mat);
-
-  if (this->non_local_manager) {
-    this->non_local_manager->initialize();
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-Int SolidMechanicsModel::getInternalIndexFromID(const ID & id) const {
-  AKANTU_DEBUG_IN();
-
-  auto it = materials.begin();
-  auto end = materials.end();
-
-  for (; it != end; ++it) {
-    if ((*it)->getID() == id) {
-      AKANTU_DEBUG_OUT();
-      return (it - materials.begin());
-    }
-  }
-
-  AKANTU_DEBUG_OUT();
-  return -1;
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::reassignMaterial() {
-  AKANTU_DEBUG_IN();
-
-  std::vector<Array<Element>> element_to_add(materials.size());
-  std::vector<Array<Element>> element_to_remove(materials.size());
-
-  for_each_element(
-      mesh,
-      [&](auto && element) {
-        auto old_material = material_index(element);
-        auto new_material = (*material_selector)(element);
-        if (old_material != new_material) {
-          element_to_add[new_material].push_back(element);
-          element_to_remove[old_material].push_back(element);
-        }
-      },
-      _spatial_dimension = spatial_dimension);
-
-  for (auto && data : enumerate(materials)) {
-    auto mat_index = std::get<0>(data);
-    auto & mat = *std::get<1>(data);
-
-    mat.removeElements(element_to_remove[mat_index]);
-    mat.addElements(element_to_add[mat_index]);
-  }
-
-  AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::applyEigenGradU(
-    const Matrix<Real> & prescribed_eigen_grad_u, const ID & material_name,
-    const GhostType ghost_type) {
-  AKANTU_DEBUG_ASSERT(prescribed_eigen_grad_u.size() ==
-                          spatial_dimension * spatial_dimension,
-                      "The prescribed grad_u is not of the good size");
-  for (auto & material : materials) {
-    if (material->getName() == material_name) {
-      material->applyEigenGradU(prescribed_eigen_grad_u, ghost_type);
-    }
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-
-} // namespace akantu
diff --git a/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh b/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh
deleted file mode 100644
index 441ed43fd..000000000
--- a/src/model/solid_mechanics/solid_mechanics_model_tmpl.hh
+++ /dev/null
@@ -1,63 +0,0 @@
-/**
- * @file   solid_mechanics_model_tmpl.hh
- *
- * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
- * @author Dana Christen <dana.christen@epfl.ch>
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
- *
- * @date creation: Fri Jun 18 2010
- * @date last modification: Fri Mar 26 2021
- *
- * @brief  template part of solid mechanics model
- *
- *
- * @section LICENSE
- *
- * Copyright (©) 2010-2021 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/>.
- *
- */
-
-/* -------------------------------------------------------------------------- */
-#include "material.hh"
-#include "solid_mechanics_model.hh"
-/* -------------------------------------------------------------------------- */
-
-#ifndef AKANTU_SOLID_MECHANICS_MODEL_TMPL_HH_
-#define AKANTU_SOLID_MECHANICS_MODEL_TMPL_HH_
-
-namespace akantu {
-
-#define FWD(...) ::std::forward<decltype(__VA_ARGS__)>(__VA_ARGS__)
-
-/* -------------------------------------------------------------------------- */
-template <typename Operation>
-void SolidMechanicsModel::splitByMaterial(const Array<Element> & elements,
-                                          Operation && op) const {
-  std::vector<Array<Element>> elements_per_mat(materials.size());
-  this->splitElementByMaterial(elements, elements_per_mat);
-
-  for (auto && mat : zip(materials, elements_per_mat)) {
-    FWD(op)(FWD(*std::get<0>(mat)), FWD(std::get<1>(mat)));
-  }
-}
-
-#undef FWD
-/* -------------------------------------------------------------------------- */
-
-} // namespace akantu
-
-#endif /* AKANTU_SOLID_MECHANICS_MODEL_TMPL_HH_ */