diff --git a/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.cpp b/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.cpp
index a19b011d6f2b4b32c5dc072e224cd5846995cfb5..a8ec6dc510e1c5db25cb48bec7684553da1ea096 100644
--- a/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.cpp
+++ b/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.cpp
@@ -62,73 +62,106 @@ private:
 void getFractureMatrixDataInMesh(
         MeshLib::Mesh const& mesh,
         std::vector<MeshLib::Element*>& vec_matrix_elements,
-        std::vector<MeshLib::Element*>& vec_fracture_elements,
-        std::vector<MeshLib::Element*>& vec_fracture_matrix_elements,
-        std::vector<MeshLib::Node*>& vec_fracture_nodes
+        std::vector<int>& vec_fracture_mat_IDs,
+        std::vector<std::vector<MeshLib::Element*>>& vec_fracture_elements,
+        std::vector<std::vector<MeshLib::Element*>>& vec_fracture_matrix_elements,
+        std::vector<std::vector<MeshLib::Node*>>& vec_fracture_nodes
         )
 {
     IsCrackTip isCrackTip(mesh);
 
     // get vectors of matrix elements and fracture elements
     vec_matrix_elements.reserve(mesh.getNumberOfElements());
+    std::vector<MeshLib::Element*> all_fracture_elements;
     for (MeshLib::Element* e : mesh.getElements())
     {
         if (e->getDimension() == mesh.getDimension())
             vec_matrix_elements.push_back(e);
         else
-            vec_fracture_elements.push_back(e);
+            all_fracture_elements.push_back(e);
     }
     DBUG("-> found total %d matrix elements and %d fracture elements",
-         vec_matrix_elements.size(), vec_fracture_elements.size());
+         vec_matrix_elements.size(), all_fracture_elements.size());
+
+    // get fracture material IDs
+    auto opt_material_ids(mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
+    if (!opt_material_ids)
+        OGS_FATAL("MaterialIDs propery vector not found in a mesh");
+    for (MeshLib::Element* e : all_fracture_elements)
+        vec_fracture_mat_IDs.push_back((*opt_material_ids)[e->getID()]);
+    std::sort(vec_fracture_mat_IDs.begin(), vec_fracture_mat_IDs.end());
+    vec_fracture_mat_IDs.erase(
+        std::unique(vec_fracture_mat_IDs.begin(), vec_fracture_mat_IDs.end()),
+        vec_fracture_mat_IDs.end());
+    DBUG("-> found %d fracture material groups", vec_fracture_mat_IDs.size());
+
+    // create a vector of fracture elements for each group
+    vec_fracture_elements.resize(vec_fracture_mat_IDs.size());
+    for (unsigned frac_id=0; frac_id<vec_fracture_mat_IDs.size(); frac_id++)
+    {
+        const auto frac_mat_id = vec_fracture_mat_IDs[frac_id];
+        std::vector<MeshLib::Element*> &vec_elements = vec_fracture_elements[frac_id];
+        for (MeshLib::Element* e : all_fracture_elements)
+        {
+            if ((*opt_material_ids)[e->getID()] == frac_mat_id)
+                vec_elements.push_back(e);
+        }
+        DBUG("-> found %d elements on the fracture %d", vec_elements.size(), frac_id);
+    }
 
     // get a vector of fracture nodes
-    for (MeshLib::Element* e : vec_fracture_elements)
+    vec_fracture_nodes.resize(vec_fracture_mat_IDs.size());
+    for (unsigned frac_id=0; frac_id<vec_fracture_mat_IDs.size(); frac_id++)
     {
-        for (unsigned i=0; i<e->getNumberOfNodes(); i++)
+        std::vector<MeshLib::Node*> &vec_nodes = vec_fracture_nodes[frac_id];
+        for (MeshLib::Element* e : vec_fracture_elements[frac_id])
         {
-            if (isCrackTip(*e->getNode(i)))
-                continue;
-            vec_fracture_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
+            for (unsigned i=0; i<e->getNumberOfNodes(); i++)
+            {
+                if (isCrackTip(*e->getNode(i)))
+                    continue;
+                vec_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
+            }
         }
+        std::sort(vec_nodes.begin(), vec_nodes.end(),
+            [](MeshLib::Node* node1, MeshLib::Node* node2) { return (node1->getID() < node2->getID()); }
+            );
+        vec_nodes.erase(std::unique(vec_nodes.begin(), vec_nodes.end()), vec_nodes.end());
+        DBUG("-> found %d nodes on the fracture %d", vec_nodes.size(), frac_id);
     }
-    std::sort(vec_fracture_nodes.begin(), vec_fracture_nodes.end(),
-        [](MeshLib::Node* node1, MeshLib::Node* node2) { return (node1->getID() < node2->getID()); }
-        );
-    vec_fracture_nodes.erase(
-                std::unique(vec_fracture_nodes.begin(), vec_fracture_nodes.end()),
-                vec_fracture_nodes.end());
-    DBUG("-> found %d nodes on the fracture", vec_fracture_nodes.size());
 
     // create a vector fracture elements and connected matrix elements,
     // which are passed to a DoF table
-    // first, collect matrix elements
-    for (MeshLib::Element *e : vec_fracture_elements)
+    for (auto fracture_elements : vec_fracture_elements)
     {
-        for (unsigned i=0; i<e->getNumberOfBaseNodes(); i++)
+        std::vector<MeshLib::Element*> vec_ele;
+        // first, collect matrix elements
+        for (MeshLib::Element*e : fracture_elements)
         {
-            MeshLib::Node const* node = e->getNode(i);
-            if (isCrackTip(*node))
-                continue;
-            for (unsigned j=0; j<node->getNumberOfElements(); j++)
+            for (unsigned i=0; i<e->getNumberOfBaseNodes(); i++)
             {
-                // only matrix elements
-                if (node->getElement(j)->getDimension() == mesh.getDimension()-1)
+                MeshLib::Node const* node = e->getNode(i);
+                if (isCrackTip(*node))
                     continue;
-                vec_fracture_matrix_elements.push_back(const_cast<MeshLib::Element*>(node->getElement(j)));
+                for (unsigned j=0; j<node->getNumberOfElements(); j++)
+                {
+                    // only matrix elements
+                    if (node->getElement(j)->getDimension() < mesh.getDimension())
+                        continue;
+                    vec_ele.push_back(const_cast<MeshLib::Element*>(node->getElement(j)));
+                }
             }
         }
+        std::sort(vec_ele.begin(), vec_ele.end(),
+            [](MeshLib::Element* p1, MeshLib::Element* p2) { return (p1->getID() < p2->getID()); }
+            );
+        vec_ele.erase(std::unique(vec_ele.begin(), vec_ele.end()), vec_ele.end());
+
+        // second, append fracture elements
+        vec_ele.insert(vec_ele.end(), fracture_elements.begin(), fracture_elements.end());
+
+        vec_fracture_matrix_elements.push_back(vec_ele);
     }
-    std::sort(vec_fracture_matrix_elements.begin(), vec_fracture_matrix_elements.end(),
-        [](MeshLib::Element* p1, MeshLib::Element* p2) { return (p1->getID() < p2->getID()); }
-        );
-    vec_fracture_matrix_elements.erase(
-                std::unique(vec_fracture_matrix_elements.begin(), vec_fracture_matrix_elements.end()),
-                vec_fracture_matrix_elements.end());
-
-    // second, append fracture elements
-    vec_fracture_matrix_elements.insert(
-                vec_fracture_matrix_elements.end(),
-                vec_fracture_elements.begin(), vec_fracture_elements.end());
 }
 
 }  // namespace SmallDeformationWithLIE
diff --git a/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.h b/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.h
index 7667c01613fdc900041b223296f4a00181fcf036..e2168211943b01f6934870cdb8ba0542cec101b7 100644
--- a/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.h
+++ b/ProcessLib/SmallDeformationWithLIE/Common/MeshUtils.h
@@ -26,16 +26,19 @@ namespace SmallDeformationWithLIE
  * @param mesh  A mesh which includes fracture elements, i.e. lower-dimensional elements.
  * It is assumed that elements forming a fracture have a distinct material ID.
  * @param vec_matrix_elements  a vector of matrix elements
- * @param vec_fracture_elements  a vector of fracture elements
- * @param vec_fracture_matrix_elements  a vector of fracture elements and matrix elements connecting to the fracture
- * @param vec_fracture_nodes  a vector of fracture element nodes
+ * @param vec_fracture_mat_IDs  fracture material IDs found in the mesh
+ * @param vec_fracture_elements  a vector of fracture elements (grouped by fracture IDs)
+ * @param vec_fracture_matrix_elements  a vector of fracture elements and matrix elements
+ * connecting to the fracture (grouped by fracture IDs)
+ * @param vec_fracture_nodes  a vector of fracture element nodes (grouped by fracture IDs)
  */
 void getFractureMatrixDataInMesh(
         MeshLib::Mesh const& mesh,
         std::vector<MeshLib::Element*>& vec_matrix_elements,
-        std::vector<MeshLib::Element*>& vec_fracture_elements,
-        std::vector<MeshLib::Element*>& vec_fracture_matrix_elements,
-        std::vector<MeshLib::Node*>& vec_fracture_nodes
+        std::vector<int>& vec_fracture_mat_IDs,
+        std::vector<std::vector<MeshLib::Element*>>& vec_fracture_elements,
+        std::vector<std::vector<MeshLib::Element*>>& vec_fracture_matrix_elements,
+        std::vector<std::vector<MeshLib::Node*>>& vec_fracture_nodes
         );
 
 }  // namespace SmallDeformationWithLIE