diff --git a/python/py_group_manager.cc b/python/py_group_manager.cc
index 02bd999da..349b757ec 100644
--- a/python/py_group_manager.cc
+++ b/python/py_group_manager.cc
@@ -1,178 +1,181 @@
 /**
  * @file   py_group_manager.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Mon Dec 02 2019
  *
  * @brief  pybind11 interface to GroupManager, ElementGroup and NodeGroup
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-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 "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <element_group.hh>
 #include <node_group.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 void register_group_manager(py::module & mod) {
   /* ------------------------------------------------------------------------ */
   py::class_<NodeGroup>(mod, "NodeGroup")
       .def(
           "getNodes",
           [](NodeGroup & self) -> decltype(auto) { return self.getNodes(); },
           py::return_value_policy::reference)
       .def("__len__", &NodeGroup::size)
       .def(
           "__iter__",
           [](const NodeGroup & self) {
             return py::make_iterator(self.begin(), self.end());
           },
           py::keep_alive<0, 1>())
       .def("__contains__",
            [](const NodeGroup & self, UInt node) {
              return self.find(node) != UInt(-1);
            })
       .def("getName", &NodeGroup::getName)
       .def("clear", &NodeGroup::clear)
       .def("empty", &NodeGroup::empty)
       .def("append", &NodeGroup::append)
       .def("add", &NodeGroup::add, py::arg("node"),
            py::arg("check_for_duplicate") = true)
       .def("remove", &NodeGroup::add);
 
   /* ------------------------------------------------------------------------ */
   py::class_<ElementGroup>(mod, "ElementGroup")
       .def(
           "getNodeGroup",
           [](ElementGroup & self) -> decltype(auto) {
             return self.getNodeGroup();
           },
           py::return_value_policy::reference)
       .def("getName", &ElementGroup::getName)
       .def(
           "getElements",
           [](ElementGroup & self) -> decltype(auto) {
             return self.getElements();
           },
           py::return_value_policy::reference)
       .def(
           "getNodeGroup",
           [](ElementGroup & self) -> decltype(auto) {
             return self.getNodeGroup();
           },
           py::return_value_policy::reference)
       .def("__len__", [](const ElementGroup & self) { return self.size(); })
       .def("clear", [](ElementGroup & self) { self.clear(); })
       .def("empty", &ElementGroup::empty)
       .def("append", &ElementGroup::append)
+      .def("optimize", &ElementGroup::optimize)
       .def(
           "add",
           [](ElementGroup & self, const Element & element, bool add_nodes,
              bool check_for_duplicate) {
             self.add(element, add_nodes, check_for_duplicate);
           },
           py::arg("element"), py::arg("add_nodes") = false,
           py::arg("check_for_duplicate") = true)
       .def("fillFromNodeGroup", &ElementGroup::fillFromNodeGroup)
       .def("addDimension", &ElementGroup::addDimension);
 
   /* ------------------------------------------------------------------------ */
   py::class_<GroupManager>(mod, "GroupManager")
       .def(
           "getElementGroup",
           [](GroupManager & self, const std::string & name) -> decltype(auto) {
             return self.getElementGroup(name);
           },
           py::return_value_policy::reference)
       .def("iterateElementGroups",
            [](GroupManager & self) -> decltype(auto) {
              std::vector<std::reference_wrapper<ElementGroup>> groups;
              for (auto & group : self.iterateElementGroups()) {
                groups.emplace_back(group);
              }
              return groups;
            })
       .def("iterateNodeGroups",
            [](GroupManager & self) -> decltype(auto) {
              std::vector<std::reference_wrapper<NodeGroup>> groups;
              for (auto & group : self.iterateNodeGroups()) {
                groups.emplace_back(group);
              }
              return groups;
            })
       .def("createNodeGroup", &GroupManager::createNodeGroup,
            py::return_value_policy::reference)
       .def(
           "createElementGroup",
           [](GroupManager & self, const std::string & id,
-             UInt spatial_dimension, bool b) -> decltype(auto) {
-            return self.createElementGroup(id, spatial_dimension, b);
+             UInt spatial_dimension, bool replace_group) -> decltype(auto) {
+            return self.createElementGroup(id, spatial_dimension,
+                                           replace_group);
           },
-          py::return_value_policy::reference)
+          py::arg("id"), py::arg("spatial_dimension"),
+          py::arg("replace_group") = false, py::return_value_policy::reference)
       .def("createGroupsFromMeshDataUInt",
            &GroupManager::createGroupsFromMeshData<UInt>)
       .def("createElementGroupFromNodeGroup",
            &GroupManager::createElementGroupFromNodeGroup, py::arg("name"),
            py::arg("node_group"), py::arg("dimension") = _all_dimensions)
       .def(
           "getNodeGroup",
           [](GroupManager & self, const std::string & name) -> decltype(auto) {
             return self.getNodeGroup(name);
           },
           py::return_value_policy::reference)
       .def(
           "nodeGroups",
           [](GroupManager & self) {
             std::vector<NodeGroup *> groups;
             for (auto & g : self.iterateNodeGroups()) {
               groups.push_back(&g);
             }
             return groups;
           },
           py::return_value_policy::reference)
       .def(
           "elementGroups",
           [](GroupManager & self) {
             std::vector<ElementGroup *> groups;
             for (auto & g : self.iterateElementGroups()) {
               groups.push_back(&g);
             }
             return groups;
           },
           py::return_value_policy::reference)
       .def("createBoundaryGroupFromGeometry",
            &GroupManager::createBoundaryGroupFromGeometry);
 }
 
 } // namespace akantu
diff --git a/python/py_mesh.cc b/python/py_mesh.cc
index c371d9c6a..33a327a62 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,237 +1,239 @@
 /**
  * @file   py_mesh.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Philip Mueller <philip.paul.mueller@bluemail.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Mon Mar 15 2021
  *
  * @brief  pybind11 interface to Mesh
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-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_config.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <mesh.hh>
 #include <mesh_accessor.hh>
 #include <mesh_utils.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 namespace {
   /* ------------------------------------------------------------------------ */
   template <typename T>
   decltype(auto) register_element_type_map_array(py::module & mod,
                                                  const std::string & name) {
     return py::class_<ElementTypeMapArray<T>,
                       std::shared_ptr<ElementTypeMapArray<T>>>(
                mod, ("ElementTypeMapArray" + name).c_str())
         .def(py::init<const ID &, const ID &>(),
              py::arg("id") = "by_element_type_array",
              py::arg("parent_id") = "no_parent")
         .def(
             "__call__",
             [](ElementTypeMapArray<T> & self, ElementType type,
                GhostType ghost_type) -> decltype(auto) {
               return self(type, ghost_type);
             },
             py::arg("type"), py::arg("ghost_type") = _not_ghost,
             py::return_value_policy::reference, py::keep_alive<0, 1>())
         .def(
             "elementTypes",
             [](ElementTypeMapArray<T> & self, UInt _dim, GhostType _ghost_type,
                ElementKind _kind) -> std::vector<ElementType> {
               auto types = self.elementTypes(_dim, _ghost_type, _kind);
               std::vector<ElementType> _types;
               for (auto && t : types) {
                 _types.push_back(t);
               }
               return _types;
             },
             py::arg("dim") = _all_dimensions,
             py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular)
         .def(
             "initialize",
             [](ElementTypeMapArray<T> & self, const Mesh & mesh,
                GhostType ghost_type = _casper, UInt nb_component = 1,
                UInt spatial_dimension = UInt(-2),
                ElementKind element_kind = _ek_not_defined,
                bool with_nb_element = false,
                bool with_nb_nodes_per_element = false, T default_value = T(),
                bool do_not_default = false) {
               self.initialize(
                   mesh, _ghost_type = ghost_type, _nb_component = nb_component,
                   _spatial_dimension = (spatial_dimension == UInt(-2)
                                             ? mesh.getSpatialDimension()
                                             : spatial_dimension),
                   _element_kind = element_kind,
                   _with_nb_element = with_nb_element,
                   _with_nb_nodes_per_element = with_nb_nodes_per_element,
                   _default_value = default_value,
                   _do_not_default = do_not_default);
             },
             py::arg("mesh"), py::arg("ghost_type") = _casper,
             py::arg("nb_component") = 1,
             py::arg("spatial_dimension") = UInt(-2),
             py::arg("element_kind") = _ek_not_defined,
             py::arg("with_nb_element") = false,
             py::arg("with_nb_nodes_per_element") = false,
             py::arg("default_value") = T(), py::arg("do_not_default") = false);
   }
 } // namespace
 /* -------------------------------------------------------------------------- */
 void register_mesh(py::module & mod) {
   py::class_<Mesh::PeriodicSlaves>(mod, "PeriodicSlaves")
       .def(
           "__iter__",
           [](Mesh::PeriodicSlaves & _this) {
             py::make_iterator(_this.begin(), _this.end());
           },
           py::keep_alive<0, 1>());
 
   py::class_<MeshData>(mod, "MeshData")
       .def(
           "getElementalDataUInt",
           [](MeshData & _this, const ID & name) -> decltype(auto) {
             return _this.getElementalData<UInt>(name);
           },
           py::return_value_policy::reference)
       .def(
           "getElementalDataReal",
           [](MeshData & _this, const ID & name) -> decltype(auto) {
             return _this.getElementalData<Real>(name);
           },
           py::return_value_policy::reference);
 
   py::class_<Mesh, GroupManager, Dumpable, MeshData>(mod, "Mesh",
                                                      py::multiple_inheritance())
       .def(py::init<UInt, const ID &>(), py::arg("spatial_dimension"),
            py::arg("id") = "mesh")
       .def("read", &Mesh::read, py::arg("filename"),
            py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file")
       .def(
           "getNodes",
           [](Mesh & self) -> decltype(auto) { return self.getNodes(); },
           py::return_value_policy::reference)
       .def("getNbNodes", &Mesh::getNbNodes)
       .def(
           "getConnectivity",
           [](Mesh & self, ElementType type) -> decltype(auto) {
             return self.getConnectivity(type);
           },
           py::return_value_policy::reference)
       .def(
           "addConnectivityType",
           [](Mesh & self, ElementType type, GhostType ghost_type) -> void {
             self.addConnectivityType(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def("distribute", [](Mesh & self) { self.distribute(); })
       .def("fillNodesToElements", &Mesh::fillNodesToElements,
            py::arg("dimension") = _all_dimensions)
       .def("getAssociatedElements",
            [](Mesh & self, const UInt & node, py::list list) {
              Array<Element> elements;
              self.getAssociatedElements(node, elements);
              for (auto && element : elements) {
                list.append(element);
              }
            })
       .def("makePeriodic",
            [](Mesh & self, const SpatialDirection & direction) {
              self.makePeriodic(direction);
            })
       .def(
           "getNbElement",
           [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type,
              ElementKind kind) {
             return self.getNbElement(spatial_dimension, ghost_type, kind);
           },
           py::arg("spatial_dimension") = _all_dimensions,
           py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined)
       .def(
           "getNbElement",
           [](Mesh & self, ElementType type, GhostType ghost_type) {
             return self.getNbElement(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def_static(
           "getSpatialDimension",
           [](ElementType & type) { return Mesh::getSpatialDimension(type); })
       .def(
           "getDataReal",
           [](Mesh & _this, const ID & name, ElementType type,
              GhostType ghost_type) -> decltype(auto) {
             return _this.getData<Real>(name, type, ghost_type);
           },
           py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "hasDataReal",
           [](Mesh & _this, const ID & name, ElementType type,
              GhostType ghost_type) -> bool {
             return _this.hasData<Real>(name, type, ghost_type);
           },
           py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); })
       .def("getPeriodicMaster", &Mesh::getPeriodicMaster)
       .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves)
       .def("isPeriodicSlave", &Mesh::isPeriodicSlave)
-      .def("isPeriodicMaster", &Mesh::isPeriodicMaster);
+      .def("isPeriodicMaster", &Mesh::isPeriodicMaster)
+      .def("initMeshFacets", &Mesh::initMeshFacets,
+           py::arg("id") = "mesh_facets", py::return_value_policy::reference);
 
   /* ------------------------------------------------------------------------ */
   py::class_<MeshUtils>(mod, "MeshUtils")
       .def_static("buildFacets", &MeshUtils::buildFacets);
 
   py::class_<MeshAccessor>(mod, "MeshAccessor")
       .def(py::init<Mesh &>(), py::arg("mesh"))
       .def(
           "resizeConnectivity",
           [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt)
               -> void { self.resizeConnectivity(new_size, type, gt); },
           py::arg("new_size"), py::arg("type"),
           py::arg("ghost_type") = _not_ghost)
       .def(
           "resizeNodes",
           [](MeshAccessor & self, UInt new_size) -> void {
             self.resizeNodes(new_size);
           },
           py::arg("new_size"))
       .def("makeReady", &MeshAccessor::makeReady);
 
   register_element_type_map_array<Real>(mod, "Real");
   register_element_type_map_array<UInt>(mod, "UInt");
   // register_element_type_map_array<std::string>(mod, "String");
 }
 } // namespace akantu
diff --git a/python/py_solid_mechanics_model_cohesive.cc b/python/py_solid_mechanics_model_cohesive.cc
index 6d4ce5f3e..ba64b18c5 100644
--- a/python/py_solid_mechanics_model_cohesive.cc
+++ b/python/py_solid_mechanics_model_cohesive.cc
@@ -1,95 +1,96 @@
 /**
  * @file   py_solid_mechanics_model_cohesive.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Tue Jul 21 2020
  * @date last modification: Tue Sep 29 2020
  *
  * @brief  pybind11 interface to SolidMechanicsModelCohesive
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-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 "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <non_linear_solver.hh>
 #include <solid_mechanics_model_cohesive.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 #define def_deprecated(func_name, mesg)                                        \
   def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); })
 
 #define def_function_nocopy(func_name)                                         \
   def(                                                                         \
       #func_name,                                                              \
       [](SolidMechanicsModel & self) -> decltype(auto) {                       \
         return self.func_name();                                               \
       },                                                                       \
       py::return_value_policy::reference)
 
 #define def_function(func_name)                                                \
   def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) {           \
     return self.func_name();                                                   \
   })
 
 void register_solid_mechanics_model_cohesive(py::module & mod) {
   py::class_<CohesiveElementInserter>(mod, "CohesiveElementInserter")
       .def("setLimit", &CohesiveElementInserter::setLimit);
 
   py::class_<SolidMechanicsModelCohesiveOptions, SolidMechanicsModelOptions>(
       mod, "SolidMechanicsModelCohesiveOptions")
       .def(py::init<AnalysisMethod, bool>(),
            py::arg("analysis_method") = _explicit_lumped_mass,
            py::arg("is_extrinsic") = false);
 
   py::class_<SolidMechanicsModelCohesive, SolidMechanicsModel>(
       mod, "SolidMechanicsModelCohesive")
       .def(py::init<Mesh &, UInt, const ID &>(), py::arg("mesh"),
            py::arg("spatial_dimension") = _all_dimensions,
            py::arg("id") = "solid_mechanics_model")
       .def(
           "initFull",
           [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method,
              bool is_extrinsic) {
             self.initFull(_analysis_method = analysis_method,
                           _is_extrinsic = is_extrinsic);
           },
-          py::arg("_analysis_method"), py::arg("_is_extrinsic") = false)
+          py::arg("_analysis_method") = _explicit_lumped_mass,
+          py::arg("_is_extrinsic") = false)
 
       .def("checkCohesiveStress",
            &SolidMechanicsModelCohesive::checkCohesiveStress)
       .def("getElementInserter",
            &SolidMechanicsModelCohesive::getElementInserter,
            py::return_value_policy::reference)
       .def("updateAutomaticInsertion",
            &SolidMechanicsModelCohesive::updateAutomaticInsertion);
 }
 
 } // namespace akantu
diff --git a/src/mesh/element_group.cc b/src/mesh/element_group.cc
index 380d2a8cd..849ddeb67 100644
--- a/src/mesh/element_group.cc
+++ b/src/mesh/element_group.cc
@@ -1,209 +1,205 @@
 /**
  * @file   element_group.cc
  *
  * @author Dana Christen <dana.christen@gmail.com>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Wed Nov 13 2013
  * @date last modification: Wed Dec 09 2020
  *
  * @brief  Stores information relevent to the notion of domain boundary and
  * surfaces.
  *
  *
  * @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_csr.hh"
 #include "dumpable.hh"
 #include "dumpable_inline_impl.hh"
 #include "group_manager.hh"
 #include "group_manager_inline_impl.hh"
 #include "mesh.hh"
 #include "mesh_utils.hh"
 #include <algorithm>
 #include <iterator>
 #include <sstream>
 
 #include "element_group.hh"
 /* -------------------------------------------------------------------------- */
 #include "dumper_iohelper_paraview.hh"
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 ElementGroup::ElementGroup(const std::string & group_name, const Mesh & mesh,
                            NodeGroup & node_group, UInt dimension,
                            const std::string & id)
     : mesh(mesh), name(group_name), elements("elements", id),
       node_group(node_group), dimension(dimension) {
   AKANTU_DEBUG_IN();
 
   this->registerDumper<DumperParaview>("paraview_" + group_name, group_name,
                                        true);
   this->addDumpFilteredMesh(mesh, elements, node_group.getNodes(),
                             _all_dimensions);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 ElementGroup::ElementGroup(const ElementGroup & /*other*/) = default;
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::clear() { elements.free(); }
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::clear(ElementType type, GhostType ghost_type) {
   if (elements.exists(type, ghost_type)) {
     elements(type, ghost_type).clear();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 bool ElementGroup::empty() const { return elements.empty(); }
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::append(const ElementGroup & other_group) {
   AKANTU_DEBUG_IN();
 
   node_group.append(other_group.node_group);
 
   /// loop on all element types in all dimensions
   for (auto ghost_type : ghost_types) {
     for (auto type : other_group.elementTypes(_ghost_type = ghost_type,
                      _element_kind = _ek_not_defined)) {
       const Array<UInt> & other_elem_list =
           other_group.elements(type, ghost_type);
       UInt nb_other_elem = other_elem_list.size();
 
       Array<UInt> * elem_list;
       UInt nb_elem = 0;
 
       /// create current type if doesn't exists, otherwise get information
       if (elements.exists(type, ghost_type)) {
         elem_list = &elements(type, ghost_type);
         nb_elem = elem_list->size();
       } else {
         elem_list = &(elements.alloc(0, 1, type, ghost_type));
       }
 
       /// append new elements to current list
       elem_list->resize(nb_elem + nb_other_elem);
       std::copy(other_elem_list.begin(), other_elem_list.end(),
                 elem_list->begin() + nb_elem);
-
-      /// remove duplicates
-      std::sort(elem_list->begin(), elem_list->end());
-      Array<UInt>::iterator<> end =
-          std::unique(elem_list->begin(), elem_list->end());
-      elem_list->resize(end - elem_list->begin());
     }
   }
 
+  this->optimize();
+
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::printself(std::ostream & stream, int indent) const {
   std::string space;
   for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) {
     ;
   }
 
   stream << space << "ElementGroup [" << std::endl;
   stream << space << " + name: " << name << std::endl;
   stream << space << " + dimension: " << dimension << std::endl;
   elements.printself(stream, indent + 1);
   node_group.printself(stream, indent + 1);
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::optimize() {
   // increasing the locality of data when iterating on the element of a group
   for (auto ghost_type : ghost_types) {
     for (auto type : elements.elementTypes(_ghost_type = ghost_type)) {
       Array<UInt> & els = elements(type, ghost_type);
       std::sort(els.begin(), els.end());
 
       Array<UInt>::iterator<> end = std::unique(els.begin(), els.end());
       els.resize(end - els.begin());
     }
   }
 
   node_group.optimize();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::fillFromNodeGroup() {
   CSR<Element> node_to_elem;
   MeshUtils::buildNode2Elements(this->mesh, node_to_elem, this->dimension);
 
   std::set<Element> seen;
 
   Array<UInt>::const_iterator<> itn = this->node_group.begin();
   Array<UInt>::const_iterator<> endn = this->node_group.end();
   for (; itn != endn; ++itn) {
     CSR<Element>::iterator ite = node_to_elem.begin(*itn);
     CSR<Element>::iterator ende = node_to_elem.end(*itn);
     for (; ite != ende; ++ite) {
       const Element & elem = *ite;
       if (this->dimension != _all_dimensions &&
           this->dimension != Mesh::getSpatialDimension(elem.type)) {
         continue;
       }
       if (seen.find(elem) != seen.end()) {
         continue;
       }
 
       UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(elem.type);
       Array<UInt>::const_iterator<Vector<UInt>> conn_it =
           this->mesh.getConnectivity(elem.type, elem.ghost_type)
               .begin(nb_nodes_per_element);
       const Vector<UInt> & conn = conn_it[elem.element];
 
       UInt count = 0;
       for (UInt n = 0; n < conn.size(); ++n) {
         count +=
             (this->node_group.getNodes().find(conn(n)) != UInt(-1) ? 1 : 0);
       }
 
       if (count == nb_nodes_per_element) {
         this->add(elem);
       }
 
       seen.insert(elem);
     }
   }
 
   this->optimize();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementGroup::addDimension(UInt dimension) {
   this->dimension = std::max(dimension, this->dimension);
 }
 
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu