diff --git a/src/Mod/Part/App/BRepMesh.cpp b/src/Mod/Part/App/BRepMesh.cpp
index a13251186f..1f1dd87246 100644
--- a/src/Mod/Part/App/BRepMesh.cpp
+++ b/src/Mod/Part/App/BRepMesh.cpp
@@ -25,6 +25,7 @@
 #include "PreCompiled.h"
 #ifndef _PreComp_
 #include <algorithm>
+#include <Precision.hxx>
 #endif
 
 #include "BRepMesh.h"
@@ -32,36 +33,218 @@
 
 using namespace Part;
 
-namespace Part {
+namespace {
 struct MeshVertex
 {
     Base::Vector3d p;
-    std::size_t i;
+    std::size_t i = 0;
 
     explicit MeshVertex(const Base::Vector3d& p)
-        : p(p), i(0)
+        : p(p)
     {
     }
 
     Base::Vector3d toPoint() const
-    { return p; }
+    {
+        return p;
+    }
 
     bool operator < (const MeshVertex &v) const
     {
-        if (fabs ( p.x - v.p.x) >= epsilon)
+        if (p.x != v.p.x) {
             return p.x < v.p.x;
-        if (fabs ( p.y - v.p.y) >= epsilon)
+        }
+        if (p.y != v.p.y) {
             return p.y < v.p.y;
-        if (fabs ( p.z - v.p.z) >= epsilon)
+        }
+        if (p.z != v.p.z) {
             return p.z < v.p.z;
-        return false; // points are considered to be equal
+        }
+
+        // points are equal
+        return false;
+    }
+};
+
+class MergeVertex
+{
+public:
+    using Facet = BRepMesh::Facet;
+
+    MergeVertex(std::vector<Base::Vector3d> points,
+                std::vector<Facet> faces,
+                double tolerance)
+        : points{std::move(points)}
+        , faces{std::move(faces)}
+        , tolerance{tolerance}
+    {
+        setDefaultMap();
+        check();
+    }
+
+    bool hasDuplicatedPoints() const
+    {
+        return duplicatedPoints > 0;
+    }
+
+    void mergeDuplicatedPoints()
+    {
+        if (!hasDuplicatedPoints()) {
+            return;
+        }
+
+        redirectPointIndex();
+        auto degreeMap = getPointDegrees();
+        decrementPointIndex(degreeMap);
+        removeUnusedPoints(degreeMap);
+        reset();
+    }
+
+    std::vector<Base::Vector3d> getPoints() const
+    {
+        return points;
+    }
+
+    std::vector<Facet> getFacets() const
+    {
+        return faces;
     }
 
 private:
-    static const double epsilon;
-};
+    void setDefaultMap()
+    {
+        // by default map point index to itself
+        mapPointIndex.resize(points.size());
+        std::generate(mapPointIndex.begin(),
+                      mapPointIndex.end(),
+                      Base::iotaGen<std::size_t>(0));
+    }
 
-const double MeshVertex::epsilon = 10 * std::numeric_limits<double>::epsilon();
+    void reset()
+    {
+        mapPointIndex.clear();
+        duplicatedPoints = 0;
+    }
+
+    void check()
+    {
+        using VertexIterator = std::vector<Base::Vector3d>::const_iterator;
+
+        double tol3d = tolerance;
+        auto vertexLess = [tol3d](const VertexIterator& v1,
+                                  const VertexIterator& v2)
+        {
+            if (fabs(v1->x - v2->x) >= tol3d) {
+                return v1->x < v2->x;
+            }
+            if (fabs(v1->y - v2->y) >= tol3d) {
+                return v1->y < v2->y;
+            }
+            if (fabs(v1->z - v2->z) >= tol3d) {
+                return v1->z < v2->z;
+            }
+            return false;  // points are considered to be equal
+        };
+        auto vertexEqual = [&](const VertexIterator& v1,
+                               const VertexIterator& v2)
+        {
+            if (vertexLess(v1, v2)) {
+                return false;
+            }
+            if (vertexLess(v2, v1)) {
+                return false;
+            }
+            return true;
+        };
+
+        std::vector<VertexIterator> vertices;
+        vertices.reserve(points.size());
+        for (auto it = points.cbegin(); it != points.cend(); ++it) {
+            vertices.push_back(it);
+        }
+
+        std::sort(vertices.begin(), vertices.end(), vertexLess);
+
+        auto next = vertices.begin();
+        while (next != vertices.end()) {
+            next = std::adjacent_find(next, vertices.end(), vertexEqual);
+            if (next != vertices.end()) {
+                auto first = next;
+                std::size_t first_index = *first - points.begin();
+                ++next;
+                while (next != vertices.end() && vertexEqual(*first, *next)) {
+                    std::size_t next_index = *next - points.begin();
+                    mapPointIndex[next_index] = first_index;
+                    ++duplicatedPoints;
+                    ++next;
+                }
+            }
+        }
+    }
+
+    void redirectPointIndex()
+    {
+        for (auto& face : faces) {
+            face.I1 = int(mapPointIndex[face.I1]);
+            face.I2 = int(mapPointIndex[face.I2]);
+            face.I3 = int(mapPointIndex[face.I3]);
+        }
+    }
+
+    std::vector<std::size_t> getPointDegrees() const
+    {
+        std::vector<std::size_t> degreeMap;
+        degreeMap.resize(points.size());
+        for (const auto& face : faces) {
+            degreeMap[face.I1]++;
+            degreeMap[face.I2]++;
+            degreeMap[face.I3]++;
+        }
+
+        return degreeMap;
+    }
+
+    void decrementPointIndex(const std::vector<std::size_t>& degreeMap)
+    {
+        std::vector<std::size_t> decrements;
+        decrements.resize(points.size());
+
+        std::size_t decr = 0;
+        for (std::size_t pos = 0; pos < points.size(); pos++) {
+            decrements[pos] = decr;
+            if (degreeMap[pos] == 0) {
+                decr++;
+            }
+        }
+
+        for (auto& face : faces) {
+            face.I1 -= int(decrements[face.I1]);
+            face.I2 -= int(decrements[face.I2]);
+            face.I3 -= int(decrements[face.I3]);
+        }
+    }
+
+    void removeUnusedPoints(const std::vector<std::size_t>& degreeMap)
+    {
+        // remove unreferenced points
+        std::vector<Base::Vector3d> new_points;
+        new_points.reserve(points.size() - duplicatedPoints);
+        for (std::size_t pos = 0; pos < points.size(); ++pos) {
+            if (degreeMap[pos] > 0) {
+                new_points.push_back(points[pos]);
+            }
+        }
+
+        points.swap(new_points);
+    }
+
+private:
+    std::vector<Base::Vector3d> points;
+    std::vector<Facet> faces;
+    double tolerance = 0.0;
+    std::size_t duplicatedPoints = 0;
+    std::vector<std::size_t> mapPointIndex;
+};
 
 }
 
@@ -115,6 +298,13 @@ void BRepMesh::getFacesFromDomains(const std::vector<Domain>& domains,
         meshPoints[vertex.i] = vertex.toPoint();
     }
     points.swap(meshPoints);
+
+    MergeVertex merge(points, faces, Precision::Confusion());
+    if (merge.hasDuplicatedPoints()) {
+        merge.mergeDuplicatedPoints();
+        points = merge.getPoints();
+        faces = merge.getFacets();
+    }
 }
 
 std::vector<BRepMesh::Segment> BRepMesh::createSegments() const
diff --git a/tests/src/Mod/Part/App/BRepMesh.cpp b/tests/src/Mod/Part/App/BRepMesh.cpp
new file mode 100644
index 0000000000..02683fbf58
--- /dev/null
+++ b/tests/src/Mod/Part/App/BRepMesh.cpp
@@ -0,0 +1,176 @@
+// SPDX-License-Identifier: LGPL-2.1-or-later
+
+#include "gtest/gtest.h"
+#include "Mod/Part/App/BRepMesh.h"
+
+// NOLINTBEGIN
+class BRepMeshTest: public ::testing::Test
+{
+protected:
+    void SetUp() override
+    {}
+
+    void TearDown() override
+    {}
+
+    std::vector<Part::BRepMesh::Domain> getNoDomains() const
+    {
+        std::vector<Part::BRepMesh::Domain> domains;
+        return domains;
+    }
+
+    std::vector<Part::BRepMesh::Domain> getEmptyDomains() const
+    {
+        Part::BRepMesh::Domain domain;
+        std::vector<Part::BRepMesh::Domain> domains;
+        domains.push_back(domain);
+        domains.push_back(domain);
+        return domains;
+    }
+
+    std::vector<Part::BRepMesh::Domain> getConnectedDomains() const
+    {
+        Part::BRepMesh::Domain domain1;
+        domain1.points.emplace_back(0, 0, 0);
+        domain1.points.emplace_back(10, 0, 0);
+        domain1.points.emplace_back(10, 10, 0);
+        domain1.points.emplace_back(0, 10, 0);
+
+        {
+            Part::BRepMesh::Facet f1;
+            f1.I1 = 0;
+            f1.I2 = 1;
+            f1.I3 = 2;
+            domain1.facets.emplace_back(f1);
+        }
+        {
+            Part::BRepMesh::Facet f2;
+            f2.I1 = 0;
+            f2.I2 = 2;
+            f2.I3 = 3;
+            domain1.facets.emplace_back(f2);
+        }
+
+        Part::BRepMesh::Domain domain2;
+        domain2.points.emplace_back(0, 0, 0);
+        domain2.points.emplace_back(0, 10, 0);
+        domain2.points.emplace_back(0, 10, 10);
+        domain2.points.emplace_back(0, 0, 10);
+
+        {
+            Part::BRepMesh::Facet f1;
+            f1.I1 = 0;
+            f1.I2 = 1;
+            f1.I3 = 2;
+            domain2.facets.emplace_back(f1);
+        }
+        {
+            Part::BRepMesh::Facet f2;
+            f2.I1 = 0;
+            f2.I2 = 2;
+            f2.I3 = 3;
+            domain2.facets.emplace_back(f2);
+        }
+
+        std::vector<Part::BRepMesh::Domain> domains;
+        domains.push_back(domain1);
+        domains.push_back(domain2);
+        return domains;
+    }
+
+    std::vector<Part::BRepMesh::Domain> getUnconnectedDomains() const
+    {
+        double eps = 1.0e-10;
+        Part::BRepMesh::Domain domain1;
+        domain1.points.emplace_back(eps, eps, eps);
+        domain1.points.emplace_back(10, 0, 0);
+        domain1.points.emplace_back(10, 10, 0);
+        domain1.points.emplace_back(eps, 10, eps);
+
+        {
+            Part::BRepMesh::Facet f1;
+            f1.I1 = 0;
+            f1.I2 = 1;
+            f1.I3 = 2;
+            domain1.facets.emplace_back(f1);
+        }
+        {
+            Part::BRepMesh::Facet f2;
+            f2.I1 = 0;
+            f2.I2 = 2;
+            f2.I3 = 3;
+            domain1.facets.emplace_back(f2);
+        }
+
+        Part::BRepMesh::Domain domain2;
+        domain2.points.emplace_back(0, 0, 0);
+        domain2.points.emplace_back(0, 10, 0);
+        domain2.points.emplace_back(0, 10, 10);
+        domain2.points.emplace_back(0, 0, 10);
+
+        {
+            Part::BRepMesh::Facet f1;
+            f1.I1 = 0;
+            f1.I2 = 1;
+            f1.I3 = 2;
+            domain2.facets.emplace_back(f1);
+        }
+        {
+            Part::BRepMesh::Facet f2;
+            f2.I1 = 0;
+            f2.I2 = 2;
+            f2.I3 = 3;
+            domain2.facets.emplace_back(f2);
+        }
+
+        std::vector<Part::BRepMesh::Domain> domains;
+        domains.push_back(domain1);
+        domains.push_back(domain2);
+        return domains;
+    }
+};
+
+TEST_F(BRepMeshTest, testNoDomains)
+{
+    std::vector<Base::Vector3d> points;
+    std::vector<Part::BRepMesh::Facet> faces;
+    Part::BRepMesh brepMesh;
+    brepMesh.getFacesFromDomains(getNoDomains(), points, faces);
+
+    EXPECT_TRUE(points.empty());
+    EXPECT_TRUE(faces.empty());
+}
+
+TEST_F(BRepMeshTest, testEmptyDomains)
+{
+    std::vector<Base::Vector3d> points;
+    std::vector<Part::BRepMesh::Facet> faces;
+    Part::BRepMesh brepMesh;
+    brepMesh.getFacesFromDomains(getEmptyDomains(), points, faces);
+
+    EXPECT_TRUE(points.empty());
+    EXPECT_TRUE(faces.empty());
+}
+
+TEST_F(BRepMeshTest, testConnectedDomains)
+{
+    std::vector<Base::Vector3d> points;
+    std::vector<Part::BRepMesh::Facet> faces;
+    Part::BRepMesh brepMesh;
+    brepMesh.getFacesFromDomains(getConnectedDomains(), points, faces);
+
+    EXPECT_EQ(points.size(), 6);
+    EXPECT_EQ(faces.size(), 4);
+}
+
+TEST_F(BRepMeshTest, testUnconnectedDomains)
+{
+    std::vector<Base::Vector3d> points;
+    std::vector<Part::BRepMesh::Facet> faces;
+    Part::BRepMesh brepMesh;
+    brepMesh.getFacesFromDomains(getUnconnectedDomains(), points, faces);
+
+    EXPECT_EQ(points.size(), 6);
+    EXPECT_EQ(faces.size(), 4);
+}
+// NOLINTEND
diff --git a/tests/src/Mod/Part/App/CMakeLists.txt b/tests/src/Mod/Part/App/CMakeLists.txt
index 7b93f834d6..af1f3f7903 100644
--- a/tests/src/Mod/Part/App/CMakeLists.txt
+++ b/tests/src/Mod/Part/App/CMakeLists.txt
@@ -2,6 +2,7 @@
 target_sources(
     Part_tests_run
         PRIVATE
+            ${CMAKE_CURRENT_SOURCE_DIR}/BRepMesh.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/FeatureChamfer.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/FeatureCompound.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/FeatureExtrusion.cpp