diff --git a/Applications/Utils/MeshGeoTools/CMakeLists.txt b/Applications/Utils/MeshGeoTools/CMakeLists.txt
index ec344b33419e60dcb6c64a5c2d3d620679236d52..a2f044b3e7833360fdb55b6f7448117f2dd21c71 100644
--- a/Applications/Utils/MeshGeoTools/CMakeLists.txt
+++ b/Applications/Utils/MeshGeoTools/CMakeLists.txt
@@ -15,3 +15,11 @@ target_link_libraries(constructMeshesFromGeometry
     ApplicationsFileIO
     ${OGS_VTK_REQUIRED_LIBS}
 )
+
+add_executable(identifySubdomains IdentifySubdomains.cpp )
+set_target_properties(identifySubdomains PROPERTIES FOLDER Utilities)
+target_link_libraries(identifySubdomains
+    MeshGeoToolsLib
+    ApplicationsFileIO
+    ${OGS_VTK_REQUIRED_LIBS}
+)
diff --git a/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp b/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8d1d3e592945112e8f7f8c6b0afe77c1fdba4b66
--- /dev/null
+++ b/Applications/Utils/MeshGeoTools/IdentifySubdomains.cpp
@@ -0,0 +1,139 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <tclap/CmdLine.h>
+
+#include "Applications/ApplicationsLib/LogogSetup.h"
+#include "BaseLib/BuildInfo.h"
+#include "MeshGeoToolsLib/IdentifySubdomainMesh.h"
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "MeshGeoToolsLib/SearchLength.h"
+#include "MeshLib/IO/readMeshFromFile.h"
+#include "MeshLib/IO/writeMeshToFile.h"
+#include "MeshLib/Mesh.h"
+
+std::vector<std::unique_ptr<MeshLib::Mesh>> readMeshes(
+    std::vector<std::string> const& filenames)
+{
+    std::vector<std::unique_ptr<MeshLib::Mesh>> meshes;
+    meshes.reserve(filenames.size());
+
+    for (auto const& filename : filenames)
+    {
+        meshes.emplace_back(MeshLib::IO::readMeshFromFile(filename));
+    }
+    if (meshes.empty())
+    {
+        OGS_FATAL("No subdomain meshes were read.");
+    }
+    return meshes;
+}
+
+int main(int argc, char* argv[])
+{
+    ApplicationsLib::LogogSetup logog_setup;
+
+    TCLAP::CmdLine cmd(
+        "Checks if the subdomain meshes are part of the bulk mesh and writes "
+        "the 'bulk_node_ids' and the 'bulk_element_ids' in each of them.\n\n"
+        "OpenGeoSys-6 software, version " +
+            BaseLib::BuildInfo::git_describe +
+            ".\n"
+            "Copyright (c) 2012-2018, OpenGeoSys Community "
+            "(http://www.opengeosys.org)",
+        ' ', BaseLib::BuildInfo::git_describe);
+
+    TCLAP::ValueArg<bool> force_overwrite_arg(
+        "f",
+        "force",
+        "Overwrites the existing subdomain meshes. (default: do not "
+        "overwrite.)",
+        false,
+        false,
+        "true/false");
+    cmd.add(force_overwrite_arg);
+
+    TCLAP::ValueArg<std::string> output_prefix_arg(
+        "o",
+        "output_prefix",
+        "Prefix the subdomain meshes' filenames with the output prefix/path.",
+        false,
+        "",
+        "path");
+    cmd.add(output_prefix_arg);
+
+    TCLAP::ValueArg<double> search_length_arg(
+        "s",
+        "searchlength",
+        "search length determining radius for the node search algorithm. "
+        "Non-negative floating point number (default 1e-16) ",
+        false,
+        1e-16,
+        "float");
+    cmd.add(search_length_arg);
+
+    TCLAP::ValueArg<std::string> bulk_mesh_arg(
+        "m", "mesh", "the file name of the bulk mesh", true, "", "mesh file");
+    cmd.add(bulk_mesh_arg);
+
+    // All the remaining arguments are used as file names for boundary/subdomain
+    // meshes.
+    TCLAP::UnlabeledMultiArg<std::string> subdomain_meshes_filenames_arg(
+        "subdomain_meshes_filenames", "mesh file names.", true,
+        "subdomain mesh file");
+    cmd.add(subdomain_meshes_filenames_arg);
+    cmd.parse(argc, argv);
+
+    //
+    // The bulk mesh.
+    //
+    std::unique_ptr<MeshLib::Mesh> bulk_mesh{
+        MeshLib::IO::readMeshFromFile(bulk_mesh_arg.getValue())};
+
+    //
+    // Read the subdomain meshes.
+    //
+    auto const subdomain_meshes =
+        readMeshes(subdomain_meshes_filenames_arg.getValue());
+
+    //
+    // Bulk mesh node searcher.
+    //
+    auto const& mesh_node_searcher =
+        MeshGeoToolsLib::MeshNodeSearcher::getMeshNodeSearcher(
+            *bulk_mesh,
+            std::make_unique<MeshGeoToolsLib::SearchLength>(
+                search_length_arg.getValue()));
+
+    //
+    // Identify the subdomains in the bulk mesh.
+    //
+    for (auto& mesh_ptr : subdomain_meshes)
+    {
+        // If force overwrite is set or the output is to different mesh than
+        // the input mesh.
+        bool const overwrite_property_vectors =
+            force_overwrite_arg.getValue() ||
+            !output_prefix_arg.getValue().empty();
+        identifySubdomainMesh(*mesh_ptr, *bulk_mesh, mesh_node_searcher,
+                              overwrite_property_vectors);
+    }
+
+    //
+    // Output after the successful subdomain mesh identification.
+    //
+    for (auto const& mesh_ptr : subdomain_meshes)
+    {
+        MeshLib::IO::writeMeshToFile(
+            *mesh_ptr,
+            output_prefix_arg.getValue() + mesh_ptr->getName() + ".vtu");
+    }
+
+    return EXIT_SUCCESS;
+}
diff --git a/Applications/Utils/Tests.cmake b/Applications/Utils/Tests.cmake
index 5f9bd044b97b2db2b387a5ba89d96669ea9249ca..d8a6f23ccb90509913c80ed160d7ef87acd57cb4 100644
--- a/Applications/Utils/Tests.cmake
+++ b/Applications/Utils/Tests.cmake
@@ -40,6 +40,34 @@ AddTest(
     expected_post_single_joint_pcs_0_ts_1_t_1.000000.vtu post_single_joint_pcs_0_ts_1_t_1.000000.vtu u u 1e-14 1e-14
 )
 
+AddTest(
+    NAME identifySubdomains_2D_Create
+    PATH MeshGeoToolsLib/IdentifySubdomains
+    EXECUTABLE identifySubdomains
+    EXECUTABLE_ARGS -m 2D_mesh.vtu -o ${Data_BINARY_DIR}/MeshGeoToolsLib/IdentifySubdomains/new_ -- 2D_mesh_top_boundary.vtu 2D_mesh_bottom_boundary.vtu
+    REQUIREMENTS NOT OGS_USE_MPI
+    TESTER vtkdiff
+    DIFF_DATA
+    2D_mesh_top.vtu new_2D_mesh_top_boundary.vtu bulk_node_ids bulk_node_ids 0 0
+    2D_mesh_top.vtu new_2D_mesh_top_boundary.vtu bulk_element_ids bulk_element_ids 0 0
+    2D_mesh_bottom.vtu new_2D_mesh_bottom_boundary.vtu bulk_node_ids bulk_node_ids 0 0
+    2D_mesh_bottom.vtu new_2D_mesh_bottom_boundary.vtu bulk_element_ids bulk_element_ids 0 0
+)
+
+AddTest(
+    NAME identifySubdomains_2D_Check
+    PATH MeshGeoToolsLib/IdentifySubdomains
+    EXECUTABLE identifySubdomains
+    EXECUTABLE_ARGS -m 2D_mesh.vtu -o ${Data_BINARY_DIR}/MeshGeoToolsLib/IdentifySubdomains/check_ -- 2D_mesh_top.vtu 2D_mesh_bottom.vtu
+    REQUIREMENTS NOT OGS_USE_MPI
+    TESTER vtkdiff
+    DIFF_DATA
+    2D_mesh_top.vtu check_2D_mesh_top.vtu bulk_node_ids bulk_node_ids 0 0
+    2D_mesh_top.vtu check_2D_mesh_top.vtu bulk_element_ids bulk_element_ids 0 0
+    2D_mesh_bottom.vtu check_2D_mesh_bottom.vtu bulk_node_ids bulk_node_ids 0 0
+    2D_mesh_bottom.vtu check_2D_mesh_bottom.vtu bulk_element_ids bulk_element_ids 0 0
+)
+
 # Mac is producing slightly different partitioning, so the results are not
 # comparable.
 AddTest(
diff --git a/Jenkinsfile b/Jenkinsfile
index 5bfab2576e0200c2b609b92656dc87243c95e57f..f9d691e8a6fc4b197fcadc4e7e04716cc8d324e9 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -99,7 +99,7 @@ pipeline {
               //       .*potential recursive class relation.*
               recordIssues tools: [[pattern: 'build/DoxygenWarnings.log',
                 tool: [$class: 'Doxygen']]],
-                unstableTotalAll: 23
+                unstableTotalAll: 24
               dir('build/docs') { stash(name: 'doxygen') }
               publishHTML(target: [allowMissing: false, alwaysLinkToLastBuild: true,
                   keepAll: true, reportDir: 'build/docs', reportFiles: 'index.html',
diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e69a202d4b8ca26db3fe3fb64bb42c4ca287641b
--- /dev/null
+++ b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
@@ -0,0 +1,225 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include <map>
+#include <vector>
+
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+
+#include "MeshNodeSearcher.h"
+
+namespace
+{
+/// Find one-to-one mapping of subdomain nodes to the bulk mesh nodes and
+/// returns the map of ids.
+std::vector<std::size_t> identifySubdomainMeshNodes(
+    MeshLib::Mesh const& subdomain_mesh,
+    MeshGeoToolsLib::MeshNodeSearcher const& mesh_node_searcher)
+{
+    // Convert nodes pointers needed for the mesh_node_searcher algorithm.
+    auto const& nodes = subdomain_mesh.getNodes();
+    std::vector<MathLib::Point3dWithID*> subdomain_points{begin(nodes),
+                                                          end(nodes)};
+
+    auto const& bulk_node_ids =
+        mesh_node_searcher.getMeshNodeIDs(subdomain_points);
+
+    if (bulk_node_ids.size() != subdomain_mesh.getNumberOfNodes())
+    {
+        OGS_FATAL(
+            "Expected to find exactly one node in the bulk mesh for each node "
+            "of the subdomain; Found %d nodes in the bulk mesh out of %d nodes "
+            "in the subdomain.",
+            bulk_node_ids.size(), subdomain_mesh.getNumberOfNodes());
+    }
+
+    return bulk_node_ids;
+}
+
+/// Find all elements which consist of all of the given node ids.
+/// \note It is recommended to include only base nodes of elements into the
+/// search \c node_ids.
+std::vector<std::size_t> findElementsInMesh(
+    MeshLib::Mesh const& mesh, std::vector<std::size_t> const& node_ids)
+{
+    //
+    // Collect all element ids for all nodes.
+    //
+    std::vector<std::size_t> common_element_ids;
+    // Every node is connected to at least one element.
+    auto const nnodes = node_ids.size();
+    common_element_ids.reserve(nnodes);
+
+    for (auto const node_id : node_ids)
+    {
+        auto const& connected_elements = mesh.getNode(node_id)->getElements();
+        std::transform(
+            begin(connected_elements), end(connected_elements),
+            back_inserter(common_element_ids),
+            [](MeshLib::Element const* const e) { return e->getID(); });
+    }
+
+    //
+    // Count how often an element is shared by all nodes.
+    //
+    std::map<std::size_t, int> element_counts;
+    for (auto const element_id : common_element_ids)
+    {
+        element_counts[element_id]++;
+    }
+
+    //
+    // Elements which are shared by as many nodes as the input nodes are the
+    // desired elements.
+    //
+    std::vector<std::size_t> element_ids;
+    for (auto const& pair : element_counts)
+    {
+        if (pair.second == static_cast<int>(nnodes))
+        {
+            element_ids.push_back(pair.first);
+        }
+    }
+
+    return element_ids;
+}
+
+/// Tries to find unique elements' ids in the bulk mesh.
+/// The subdomain elements' nodes are mapped to the bulk mesh and are then used
+/// to identify the bulk mesh elements.
+std::vector<std::size_t> identifySubdomainMeshElements(
+    MeshLib::Mesh const& subdomain_mesh, MeshLib::Mesh const& bulk_mesh)
+{
+    auto& properties = subdomain_mesh.getProperties();
+    auto const& bulk_node_ids =
+        *properties.getPropertyVector<std::size_t>("bulk_node_ids");
+
+    // Initialize with very large value
+    std::vector<std::size_t> bulk_element_ids_map(
+        subdomain_mesh.getNumberOfElements(), -1);
+
+    for (auto* const e : subdomain_mesh.getElements())
+    {
+        std::vector<std::size_t> element_node_ids(e->getNumberOfBaseNodes());
+        for (unsigned n = 0; n < e->getNumberOfBaseNodes(); ++n)
+        {
+            element_node_ids[n] = e->getNodeIndex(n);
+        }
+        std::vector<std::size_t> element_node_ids_bulk(
+            e->getNumberOfBaseNodes());
+        std::transform(begin(element_node_ids), end(element_node_ids),
+                       begin(element_node_ids_bulk),
+                       [&bulk_node_ids](std::size_t const id) {
+                           return bulk_node_ids[id];
+                       });
+        auto const& bulk_element_ids =
+            findElementsInMesh(bulk_mesh, element_node_ids_bulk);
+
+        if (bulk_element_ids.empty())
+        {
+            ERR("No element could be found for the subdomain element %d. "
+                "Corresponding bulk mesh node ids are:",
+                e->getID());
+            for (auto const i : element_node_ids_bulk)
+                ERR("\t%d", i);
+            OGS_FATAL(
+                "Expect exactly one element to be found in the bulk mesh.");
+        }
+        if (bulk_element_ids.size() != 1)
+        {
+            ERR("%d elements found for the subdomain element %d. "
+                "Corresponding bulk mesh node ids are:",
+                bulk_element_ids.size(), e->getID());
+            for (auto const i : element_node_ids_bulk)
+                ERR("\t%d", i);
+            ERR("The found bulk mesh element ids are:");
+            for (auto const i : bulk_element_ids)
+                ERR("\t%d", i);
+            OGS_FATAL(
+                "Expect exactly one element to be found in the bulk mesh.");
+        }
+
+        bulk_element_ids_map[e->getID()] = bulk_element_ids.front();
+    }
+
+    return bulk_element_ids_map;
+}
+
+/// Updates or checks the existing mesh's property with the given values.
+void updateOrCheckExistingSubdomainProperty(
+    MeshLib::Mesh& mesh, std::string const& property_name,
+    std::vector<std::size_t> values, MeshLib::MeshItemType const mesh_item_type,
+    bool const force_overwrite)
+{
+    auto& properties = mesh.getProperties();
+    if (!properties.existsPropertyVector<std::size_t>(property_name))
+    {
+        addPropertyToMesh(mesh, property_name, mesh_item_type, 1, values);
+        return;
+    }
+
+    //
+    // Check the existing property against new values.
+    //
+    auto& original_property =
+        *properties.getPropertyVector<std::size_t>(property_name);
+    if (std::equal(begin(original_property), end(original_property),
+                   begin(values), end(values)))
+    {
+        INFO(
+            "There is already a '%s' property present in the subdomain mesh "
+            "'%s' and it is equal to the newly computed values.",
+            property_name.c_str(), mesh.getName().c_str());
+        return;
+    }
+
+    //
+    // Property differs. Notify and update if forced.
+    //
+    WARN(
+        "There is already a '%s' property present in the subdomain mesh '%s' "
+        "and it is not equal to the newly computed values.",
+        property_name.c_str(),
+        mesh.getName().c_str());
+
+    if (!force_overwrite)
+    {
+        OGS_FATAL("The force overwrite flag was not specified, exiting.");
+    }
+
+    INFO("Overwriting '%s' property.", property_name.c_str());
+    original_property.resize(values.size());
+    std::copy(begin(values), end(values), begin(original_property));
+}
+}  // namespace
+
+namespace MeshGeoToolsLib
+{
+void identifySubdomainMesh(MeshLib::Mesh& subdomain_mesh,
+                           MeshLib::Mesh const& bulk_mesh,
+                           MeshNodeSearcher const& mesh_node_searcher,
+                           bool const force_overwrite = false)
+{
+    auto const& bulk_node_ids =
+        identifySubdomainMeshNodes(subdomain_mesh, mesh_node_searcher);
+
+    updateOrCheckExistingSubdomainProperty(
+        subdomain_mesh, "bulk_node_ids", bulk_node_ids,
+        MeshLib::MeshItemType::Node, force_overwrite);
+
+    auto const& bulk_element_ids =
+        identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
+
+    updateOrCheckExistingSubdomainProperty(
+        subdomain_mesh, "bulk_element_ids", bulk_element_ids,
+        MeshLib::MeshItemType::Cell, force_overwrite);
+}
+}  // namespace MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.h b/MeshGeoToolsLib/IdentifySubdomainMesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c12013997a39383416dc63b67dbc84c6fbf06e8
--- /dev/null
+++ b/MeshGeoToolsLib/IdentifySubdomainMesh.h
@@ -0,0 +1,31 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+namespace MeshGeoToolsLib
+{
+class SearchLength;
+class MeshNodeSearcher;
+}
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace MeshGeoToolsLib
+{
+/// Geometrically finds nodes and elements of the subdomain mesh in the bulk
+/// mesh, and updates or verifies the corresponding bulk_node_ids and
+/// bulk_element_ids properties.
+///
+/// In case of unexpected results OGS_FATAL is called.
+void identifySubdomainMesh(MeshLib::Mesh& subdomain_mesh,
+                           MeshLib::Mesh const& bulk_mesh,
+                           MeshNodeSearcher const& mesh_node_searcher,
+                           bool const force_overwrite = false);
+}  // namespace MeshGeoToolsLib
diff --git a/MeshGeoToolsLib/MeshNodeSearcher.cpp b/MeshGeoToolsLib/MeshNodeSearcher.cpp
index fc19986e196685be6153ee6a05ef6cee49144eb6..7d7e5f551d8fc393374216bf771f543ba58d2de1 100644
--- a/MeshGeoToolsLib/MeshNodeSearcher.cpp
+++ b/MeshGeoToolsLib/MeshNodeSearcher.cpp
@@ -76,6 +76,41 @@ std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(
     return vec_nodes;
 }
 
+std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(
+    std::vector<MathLib::Point3dWithID*> const& points) const
+{
+    double const epsilon_radius = _search_length_algorithm->getSearchLength();
+
+    std::vector<std::size_t> node_ids;
+    node_ids.reserve(points.size());
+
+    for (auto const* const p_ptr : points)
+    {
+        auto const& p = *p_ptr;
+        std::vector<std::size_t> const ids =
+            _mesh_grid.getPointsInEpsilonEnvironment(p, epsilon_radius);
+        if (ids.empty())
+        {
+            OGS_FATAL(
+                "No nodes could be found in the mesh for point %d : (%g, %g, "
+                "%g) in %g epsilon radius in the mesh '%s'",
+                p.getID(), p[0], p[1], p[2], epsilon_radius,
+                _mesh.getName().c_str());
+        }
+        if (ids.size() != 1)
+        {
+            OGS_FATAL(
+                "Found %d nodes in the mesh for point %d : (%g, %g, %g) in %g "
+                "epsilon radius in the mesh '%s'. Expected to find exactly one "
+                "node.",
+                ids.size(), p.getID(), p[0], p[1], p[2], epsilon_radius,
+                _mesh.getName().c_str());
+        }
+        node_ids.push_back(ids.front());
+    }
+    return node_ids;
+}
+
 std::vector<std::size_t> const& MeshNodeSearcher::getMeshNodeIDsForPoint(
     GeoLib::Point const& pnt) const
 {
diff --git a/MeshGeoToolsLib/MeshNodeSearcher.h b/MeshGeoToolsLib/MeshNodeSearcher.h
index df70806611b6ab7f298ee66df7a5c921f0757f19..26d3ddc9e97e930c2c1b5720a5eef4a451589031 100644
--- a/MeshGeoToolsLib/MeshNodeSearcher.h
+++ b/MeshGeoToolsLib/MeshNodeSearcher.h
@@ -76,6 +76,14 @@ public:
     std::vector<std::size_t> getMeshNodeIDs(
         GeoLib::GeoObject const& geoObj) const;
 
+    /**
+     * Finds unique mesh nodes of each of the input points.
+     *
+     * \return a vector of mesh node ids.
+     */
+    std::vector<std::size_t> getMeshNodeIDs(
+        std::vector<MathLib::Point3dWithID*> const& points) const;
+
     /**
      * Searches for the node nearest by the given point. If there are two nodes
      * with the same distance the id of the one that was first found will be
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8f6eebbebe00a8afee59c37c410d2cdc1ec05c3d
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:578270062e2b4a11a68daa1a33434cc49aa73eedcb9debb439f66dcabce433bc
+size 55463
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_bottom.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_bottom.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..606226937dda317208c598998989f71a9620b116
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_bottom.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:823cb86a9e81ef2dbfc4375992851ca12e7ab919799ff1025916716d72105d02
+size 4700
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_bottom_boundary.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_bottom_boundary.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a40f85ea56554093b00793c513e5faab3e1c1ffd
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_bottom_boundary.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5f19b22eac51c8642cd257e27fb10476cde839e32a9c838b88da281fd8f7c6f4
+size 3561
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_top.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_top.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..69e58fc9b2b966c9c95f8f3dd92c8a0610b37f7f
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_top.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:608d488236368b497614495a69620900c8406afb63c10dd238a1cb59d40e8bb8
+size 4721
diff --git a/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_top_boundary.vtu b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_top_boundary.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..bd18fe9550b6aa47d39f9b4b3032a4655f443559
--- /dev/null
+++ b/Tests/Data/MeshGeoToolsLib/IdentifySubdomains/2D_mesh_top_boundary.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:466108735f4b863f26d47f6ac89d98602296fbf17c60f8930bb1dd3a19c61747
+size 3572
diff --git a/web/content/docs/tools/model-preparation/identifySubdomains/disc_with_hole_and_bondary.png b/web/content/docs/tools/model-preparation/identifySubdomains/disc_with_hole_and_bondary.png
new file mode 100644
index 0000000000000000000000000000000000000000..6ea6ba3ef134f1297056c389369492bd3f460509
--- /dev/null
+++ b/web/content/docs/tools/model-preparation/identifySubdomains/disc_with_hole_and_bondary.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6e4cbb161c5379f312f0f7549d792fdcb507999b8c02ccfbc1a5d447b137d135
+size 94116
diff --git a/web/content/docs/tools/model-preparation/identifySubdomains/index.pandoc b/web/content/docs/tools/model-preparation/identifySubdomains/index.pandoc
new file mode 100644
index 0000000000000000000000000000000000000000..16a43ab6e2212d68eb39868b677ff2883e26caf7
--- /dev/null
+++ b/web/content/docs/tools/model-preparation/identifySubdomains/index.pandoc
@@ -0,0 +1,51 @@
++++
+date = "2018-10-02T15:56:57+01:00"
+title = "Identitfy subdomains in bulk mesh."
+author = "Dmitri Naumov"
+
+[menu]
+  [menu.tools]
+    parent = "model-preparation"
++++
+
+## General
+
+`identifySubdomains` is a creation and check tool for the `bulk_node_ids` and
+`bulk_element_ids` mesh properties. These properties are needed for the boundary
+conditions and source terms application, defined on the subdomains of a "bulk"
+mesh.
+
+## Example
+
+Given a "bulk" mesh (Tests/Data/Mechanics/Linear/disc_with_hole.vtu) and a
+[quater cirle mesh](../quater_circle.vtu) extracted manually we want to use the
+quater circle mesh for heterogeneous boundary condition. OGS requires two
+mappings into the "bulk" mesh, one for the nodes and one for the elements.
+
+<center>
+![The figure shows a part of the "bulk" mesh with boundary element numbers, and
+the quater circle mesh shown as white line with green
+points.](disc_with_hole_and_bondary.png){width=50%}
+</center>
+
+To create this mappings we run
+```
+identifySubdomains -m Tests/Data/Mechanics/Linear/disc_with_hole.vtu -s 1e-6 -o
+new_ -- quater_circle.vtu
+```
+The tool will first try to find all unique nodes in the "bulk" mesh using search
+radius 1e-6, and create the `bulk_node_ids` mapping upon success. Then the
+`bulk_element_ids` mapping is created by finding a unique element containing all
+the nodes of the subdomain element. The output file
+[`new_quater_circle.vtu`](../new_quater_cirle.vtu) will now contain both
+mappings and is prepared for usage as a boundary condition mesh.
+
+## Notes
+
+ - The double dash is separating the subdomain meshes, so the input can have
+   multiple subdomain inputs:
+   ```
+   identifySubdomains -m bulk.vtu -- bc1.vtu bc2.vtu source_term_meshes*.vtu
+   ```
+ - The output prefix `-o` can contain a path too and is relative to the current
+   working directory.
diff --git a/web/content/docs/tools/model-preparation/identifySubdomains/new_quater_circle.vtu b/web/content/docs/tools/model-preparation/identifySubdomains/new_quater_circle.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..64653593917983705540431a28140bd755149678
--- /dev/null
+++ b/web/content/docs/tools/model-preparation/identifySubdomains/new_quater_circle.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8f1ba00fcaaa465639aec4b3931332d23e00ba18da68def9513f3f69ec8725f9
+size 3865
diff --git a/web/content/docs/tools/model-preparation/identifySubdomains/quater_circle.vtu b/web/content/docs/tools/model-preparation/identifySubdomains/quater_circle.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..64653593917983705540431a28140bd755149678
--- /dev/null
+++ b/web/content/docs/tools/model-preparation/identifySubdomains/quater_circle.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8f1ba00fcaaa465639aec4b3931332d23e00ba18da68def9513f3f69ec8725f9
+size 3865