diff --git a/source/blender/bmesh/intern/bmesh_core.cc b/source/blender/bmesh/intern/bmesh_core.cc
index bcc21562652..05d9339db8e 100644
--- a/source/blender/bmesh/intern/bmesh_core.cc
+++ b/source/blender/bmesh/intern/bmesh_core.cc
@@ -1156,7 +1156,7 @@ static bool bm_vert_is_manifold_flagged(BMVert *v, const char api_flag)
 
 /* Mid-level Topology Manipulation Functions */
 
-BMFace *BM_faces_join(BMesh *bm, BMFace **faces, int totface, const bool do_del)
+BMFace *BM_faces_join(BMesh *bm, BMFace **faces, int totface, const bool do_del, BMFace **r_double)
 {
   BMFace *f, *f_new;
 #ifdef USE_BMESH_HOLES
@@ -1168,6 +1168,12 @@ BMFace *BM_faces_join(BMesh *bm, BMFace **faces, int totface, const bool do_del)
   BMVert *v1 = nullptr, *v2 = nullptr;
   int i;
   const int cd_loop_mdisp_offset = CustomData_get_offset(&bm->ldata, CD_MDISPS);
+  BMFace *f_existing;
+
+  /* Initialize the return value if provided. This ensures it will be nullptr if the join fails. */
+  if (r_double) {
+    *r_double = nullptr;
+  }
 
   if (UNLIKELY(!totface)) {
     BMESH_ASSERT(0);
@@ -1256,6 +1262,22 @@ BMFace *BM_faces_join(BMesh *bm, BMFace **faces, int totface, const bool do_del)
     goto error;
   }
 
+  /* If a new face was created, check whether it is a double of an existing face. */
+  f_existing = BM_face_find_double(f_new);
+  if (f_existing) {
+
+    /* Return the double to the calling function if that was requested. */
+    if (r_double) {
+      *r_double = f_existing;
+    }
+
+    /* Otherwise, automatically reuse the existing face. */
+    else {
+      BM_face_kill(bm, f_new);
+      f_new = f_existing;
+    }
+  }
+
   /* copy over loop data */
   l_iter = l_first = BM_FACE_FIRST_LOOP(f_new);
   do {
diff --git a/source/blender/bmesh/intern/bmesh_core.hh b/source/blender/bmesh/intern/bmesh_core.hh
index 83e9712d5bc..b90249c3769 100644
--- a/source/blender/bmesh/intern/bmesh_core.hh
+++ b/source/blender/bmesh/intern/bmesh_core.hh
@@ -169,10 +169,35 @@ void bmesh_face_swap_data(BMFace *f_a, BMFace *f_b);
  * \note If a pair of faces share multiple edges,
  * the pair of faces will be joined at every edge.
  *
+ * \param bm: The bmesh.
+ * \param faces: An array of faces to join.
+ * \param totfaces: The length of the face array to join.
+ * \param do_del if true, remove the original faces, internal edges, and internal verts such that
+ *    they are replaced by the new face.
+ * \param r_double: A pointer to a BMFace* that is controls processing of doubled faces.
+ * - When `r_double` is nullptr:
+ *   - If a new face would be made which would double an existing face, then instead of creating a
+ *     new face, the existing face will be reused and returned instead.
+ *   - The calling function must not make ANY assumption about whether the returned BMFace* is
+ *     new, or a reused face that may already have set header flags, contain custom data, etc.
+ * - When `r_double` is a pointer to a BMFace*:
+ *   - If the new join face is not a double of an existing face, then `r_double` is set to nullptr.
+ *   - If the new join face doubles an existing face, then `r_double` is set to the existing face,
+ *     and the return value is the newly created face. The double will NOT be removed, meaning the
+ *     BMesh is in an invalid state, and the calling function must fix that inconsistency.
+ *   - If an error occurs and nullptr is returned, `r_double` will be set to nullptr as well.
+ *
  * \note this is a generic, flexible join faces function,
  * almost everything uses this, including #BM_faces_join_pair
+ *
+ * \note On callers asserting when `*r_double != nullptr`.
+ * For some callers the existing algorithm does not check for or handle double faces.
+ * This can result in invalid meshes being returned.
+ * The returned value in `r_double` should be examined and if found,
+ * the algorithm should be adjusted. Until this is changed, at least warn.
+ * This comment can be removed when all callers handle this case.
  */
-BMFace *BM_faces_join(BMesh *bm, BMFace **faces, int totface, bool do_del);
+BMFace *BM_faces_join(BMesh *bm, BMFace **faces, int totface, bool do_del, BMFace **r_double);
 /**
  * High level function which wraps both #bmesh_kernel_vert_separate and #bmesh_kernel_edge_separate
  */
diff --git a/source/blender/bmesh/intern/bmesh_mods.cc b/source/blender/bmesh/intern/bmesh_mods.cc
index 37028a765f8..b51e4cc610b 100644
--- a/source/blender/bmesh/intern/bmesh_mods.cc
+++ b/source/blender/bmesh/intern/bmesh_mods.cc
@@ -90,9 +90,16 @@ bool BM_disk_dissolve(BMesh *bm, BMVert *v)
       return false;
     }
 #else
-    if (UNLIKELY(!BM_faces_join_pair(bm, e->l, e->l->radial_next, true))) {
+    BMFace *f_double;
+
+    if (UNLIKELY(!BM_faces_join_pair(bm, e->l, e->l->radial_next, true, &f_double))) {
       return false;
     }
+
+    /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+    BLI_assert_msg(f_double == nullptr,
+                   "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
     if (UNLIKELY(!BM_vert_collapse_faces(bm, v->e, v, 1.0, true, false, true, true))) {
       return false;
     }
@@ -109,9 +116,15 @@ bool BM_disk_dissolve(BMesh *bm, BMVert *v)
 
     /* handle two-valence */
     if (e->l != e->l->radial_next) {
-      if (!BM_faces_join_pair(bm, e->l, e->l->radial_next, true)) {
+      BMFace *f_double;
+
+      if (!BM_faces_join_pair(bm, e->l, e->l->radial_next, true, &f_double)) {
         return false;
       }
+
+      /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+      BLI_assert_msg(f_double == nullptr,
+                     "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
     }
 
     return true;
@@ -126,13 +139,19 @@ bool BM_disk_dissolve(BMesh *bm, BMVert *v)
       do {
         BMFace *f = nullptr;
         if (BM_edge_is_manifold(e) && (e != baseedge) && (e != keepedge)) {
-          f = BM_faces_join_pair(bm, e->l, e->l->radial_next, true);
+          BMFace *f_double;
+
+          f = BM_faces_join_pair(bm, e->l, e->l->radial_next, true, &f_double);
           /* return if couldn't join faces in manifold
            * conditions */
           /* !disabled for testing why bad things happen */
           if (!f) {
             return false;
           }
+
+          /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+          BLI_assert_msg(f_double == nullptr,
+                         "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
         }
 
         if (f) {
@@ -154,10 +173,16 @@ bool BM_disk_dissolve(BMesh *bm, BMVert *v)
     if (e->l) {
       /* get remaining two faces */
       if (e->l != e->l->radial_next) {
+        BMFace *f_double;
+
         /* join two remaining faces */
-        if (!BM_faces_join_pair(bm, e->l, e->l->radial_next, true)) {
+        if (!BM_faces_join_pair(bm, e->l, e->l->radial_next, true, &f_double)) {
           return false;
         }
+
+        /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+        BLI_assert_msg(f_double == nullptr,
+                       "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
       }
     }
   }
@@ -165,7 +190,8 @@ bool BM_disk_dissolve(BMesh *bm, BMVert *v)
   return true;
 }
 
-BMFace *BM_faces_join_pair(BMesh *bm, BMLoop *l_a, BMLoop *l_b, const bool do_del)
+BMFace *BM_faces_join_pair(
+    BMesh *bm, BMLoop *l_a, BMLoop *l_b, const bool do_del, BMFace **r_double)
 {
   BLI_assert((l_a != l_b) && (l_a->e == l_b->e));
 
@@ -175,7 +201,7 @@ BMFace *BM_faces_join_pair(BMesh *bm, BMLoop *l_a, BMLoop *l_b, const bool do_de
   }
 
   BMFace *faces[2] = {l_a->f, l_b->f};
-  return BM_faces_join(bm, faces, 2, do_del);
+  return BM_faces_join(bm, faces, 2, do_del, r_double);
 }
 
 BMFace *BM_face_split(BMesh *bm,
@@ -369,7 +395,14 @@ BMEdge *BM_vert_collapse_faces(BMesh *bm,
     }
 
     if (faces.size() >= 2) {
-      BMFace *f2 = BM_faces_join(bm, faces.data(), faces.size(), true);
+      BMFace *f_double;
+
+      BMFace *f2 = BM_faces_join(bm, faces.data(), faces.size(), true, &f_double);
+
+      /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+      BLI_assert_msg(f_double == nullptr,
+                     "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
       if (f2) {
         BMLoop *l_a, *l_b;
 
@@ -788,9 +821,15 @@ BMEdge *BM_edge_rotate(BMesh *bm, BMEdge *e, const bool ccw, const short check_f
 
   const bool is_flipped = !BM_edge_is_contiguous(e);
 
+  BMFace *f_double;
+
   /* don't delete the edge, manually remove the edge after so we can copy its attributes */
   f = BM_faces_join_pair(
-      bm, BM_face_edge_share_loop(l1->f, e), BM_face_edge_share_loop(l2->f, e), true);
+      bm, BM_face_edge_share_loop(l1->f, e), BM_face_edge_share_loop(l2->f, e), true, &f_double);
+
+  /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+  BLI_assert_msg(f_double == nullptr,
+                 "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
 
   if (f == nullptr) {
     return nullptr;
diff --git a/source/blender/bmesh/intern/bmesh_mods.hh b/source/blender/bmesh/intern/bmesh_mods.hh
index 2d540520387..4d94a78b4a7 100644
--- a/source/blender/bmesh/intern/bmesh_mods.hh
+++ b/source/blender/bmesh/intern/bmesh_mods.hh
@@ -53,9 +53,17 @@ bool BM_disk_dissolve(BMesh *bm, BMVert *v);
  * If the windings do not match the winding of the new face will follow
  * \a l_a's winding (i.e. \a l_b will be reversed before the join).
  *
+ * \param bm: The bmesh.
+ * \param l_a: First loop of an adjacent face pair that will be joined.
+ * \param l_b: Second loop of an adjacent face pair that will be joined.
+ * \param do_del If true, remove the original faces, internal edges,
+ * and internal verts such that they are replaced by the new face.
+ * \param r_double: A pointer to a BMFace* that controls processing of doubled faces.
+ * See #BM_faces_join_pair `r_double` argument for details.
+ *
  * \return The combined face or NULL on failure.
  */
-BMFace *BM_faces_join_pair(BMesh *bm, BMLoop *l_a, BMLoop *l_b, bool do_del);
+BMFace *BM_faces_join_pair(BMesh *bm, BMLoop *l_a, BMLoop *l_b, bool do_del, BMFace **r_double);
 
 /** see: bmesh_polygon_edgenet.hh for #BM_face_split_edgenet */
 
diff --git a/source/blender/bmesh/operators/bmo_connect_concave.cc b/source/blender/bmesh/operators/bmo_connect_concave.cc
index 80b81ad098d..46a28df3935 100644
--- a/source/blender/bmesh/operators/bmo_connect_concave.cc
+++ b/source/blender/bmesh/operators/bmo_connect_concave.cc
@@ -140,8 +140,13 @@ static bool bm_face_split_by_concave(BMesh *bm,
         }
 
         if (ok) {
+          BMFace *f_double;
           BMFace *f_new, *f_pair[2] = {l_pair[0]->f, l_pair[1]->f};
-          f_new = BM_faces_join(bm, f_pair, 2, true);
+          f_new = BM_faces_join(bm, f_pair, 2, true, &f_double);
+          /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+          BLI_assert_msg(f_double == nullptr,
+                         "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
           if (f_new) {
             BMO_face_flag_enable(bm, f_new, FACE_OUT);
           }
diff --git a/source/blender/bmesh/operators/bmo_dissolve.cc b/source/blender/bmesh/operators/bmo_dissolve.cc
index f07e2eb21dc..7a4b0936946 100644
--- a/source/blender/bmesh/operators/bmo_dissolve.cc
+++ b/source/blender/bmesh/operators/bmo_dissolve.cc
@@ -189,7 +189,14 @@ void bmo_dissolve_faces_exec(BMesh *bm, BMOperator *op)
   for (Vector<BMFace *> &faces : regions) {
     const int64_t faces_len = faces.size();
 
-    BMFace *f_new = BM_faces_join(bm, faces.data(), faces_len, true);
+    BMFace *f_double;
+
+    BMFace *f_new = BM_faces_join(bm, faces.data(), faces_len, true, &f_double);
+
+    /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+    BLI_assert_msg(f_double == nullptr,
+                   "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
     if (f_new != nullptr) {
       /* Maintain the active face. */
       if (act_face && bm->act_face == nullptr) {
@@ -305,8 +312,12 @@ void bmo_dissolve_edges_exec(BMesh *bm, BMOperator *op)
     if (BM_edge_loop_pair(e, &l_a, &l_b)) {
       BMFace *f_new;
 
+      BMFace *f_double;
+
       /* join faces */
-      f_new = BM_faces_join_pair(bm, l_a, l_b, false);
+      f_new = BM_faces_join_pair(bm, l_a, l_b, false, &f_double);
+
+      /* This algorithm actually does check for double faces. */
       if (f_new && BM_face_find_double(f_new)) {
         BM_face_kill(bm, f_new);
         f_new = nullptr;
@@ -426,8 +437,12 @@ void bmo_dissolve_verts_exec(BMesh *bm, BMOperator *op)
         if (BM_edge_loop_pair(e, &l_a, &l_b)) {
           BMFace *f_new;
 
+          BMFace *f_double;
+
           /* join faces */
-          f_new = BM_faces_join_pair(bm, l_a, l_b, false);
+          f_new = BM_faces_join_pair(bm, l_a, l_b, false, &f_double);
+
+          /* This algorithm actually does check for double faces. */
           if (f_new && BM_face_find_double(f_new)) {
             BM_face_kill(bm, f_new);
             f_new = nullptr;
diff --git a/source/blender/bmesh/operators/bmo_join_triangles.cc b/source/blender/bmesh/operators/bmo_join_triangles.cc
index 2272e93ca7e..b683fb948fe 100644
--- a/source/blender/bmesh/operators/bmo_join_triangles.cc
+++ b/source/blender/bmesh/operators/bmo_join_triangles.cc
@@ -952,8 +952,15 @@ static BMFace *bm_faces_join_pair_by_edge(BMesh *bm,
   }
 #endif
 
+  BMFace *f_double;
+
   /* Join the edge and identify the face. */
-  return BM_faces_join_pair(bm, l_a, l_b, true);
+  BMFace *f = BM_faces_join_pair(bm, l_a, l_b, true, &f_double);
+  /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+  BLI_assert_msg(f_double == nullptr,
+                 "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
+  return f;
 }
 
 /** Given a mesh, convert triangles to quads. */
diff --git a/source/blender/bmesh/operators/bmo_triangulate.cc b/source/blender/bmesh/operators/bmo_triangulate.cc
index bc95bb4507d..209e2924031 100644
--- a/source/blender/bmesh/operators/bmo_triangulate.cc
+++ b/source/blender/bmesh/operators/bmo_triangulate.cc
@@ -242,7 +242,12 @@ void bmo_triangle_fill_exec(BMesh *bm, BMOperator *op)
       if (BMO_edge_flag_test(bm, e, ELE_NEW)) {
         /* in rare cases the edges face will have already been removed from the edge */
         if (LIKELY(BM_edge_is_manifold(e))) {
-          BMFace *f_new = BM_faces_join_pair(bm, e->l, e->l->radial_next, false);
+          BMFace *f_double;
+          BMFace *f_new = BM_faces_join_pair(bm, e->l, e->l->radial_next, false, &f_double);
+          /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+          BLI_assert_msg(f_double == nullptr,
+                         "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
           if (f_new) {
             BMO_face_flag_enable(bm, f_new, ELE_NEW);
             BM_edge_kill(bm, e);
diff --git a/source/blender/bmesh/tools/bmesh_decimate_collapse.cc b/source/blender/bmesh/tools/bmesh_decimate_collapse.cc
index fda4969f234..3f117006166 100644
--- a/source/blender/bmesh/tools/bmesh_decimate_collapse.cc
+++ b/source/blender/bmesh/tools/bmesh_decimate_collapse.cc
@@ -662,8 +662,14 @@ static void bm_decim_triangulate_end(BMesh *bm, const int edges_tri_tot)
     BMLoop *l_a, *l_b;
     e = edges_tri[i];
     if (BM_edge_loop_pair(e, &l_a, &l_b)) {
+      BMFace *f_double;
+
       BMFace *f_array[2] = {l_a->f, l_b->f};
-      BM_faces_join(bm, f_array, 2, false);
+      BM_faces_join(bm, f_array, 2, false, &f_double);
+      /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+      BLI_assert_msg(f_double == nullptr,
+                     "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
       if (e->l == nullptr) {
         BM_edge_kill(bm, e);
       }
diff --git a/source/blender/bmesh/tools/bmesh_decimate_dissolve.cc b/source/blender/bmesh/tools/bmesh_decimate_dissolve.cc
index 0daf9788cf5..ac3049cdbfa 100644
--- a/source/blender/bmesh/tools/bmesh_decimate_dissolve.cc
+++ b/source/blender/bmesh/tools/bmesh_decimate_dissolve.cc
@@ -348,8 +348,11 @@ void BM_mesh_decimate_dissolve_ex(BMesh *bm,
       i = BM_elem_index_get(e);
 
       if (BM_edge_is_manifold(e)) {
-        f_new = BM_faces_join_pair(bm, e->l, e->l->radial_next, false);
-
+        BMFace *f_double;
+        f_new = BM_faces_join_pair(bm, e->l, e->l->radial_next, false, &f_double);
+        /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+        BLI_assert_msg(f_double == nullptr,
+                       "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
         if (f_new) {
           BMLoop *l_first, *l_iter;
 
diff --git a/source/blender/editors/mesh/editmesh_intersect.cc b/source/blender/editors/mesh/editmesh_intersect.cc
index 8679c42bd2a..4c52110a2d6 100644
--- a/source/blender/editors/mesh/editmesh_intersect.cc
+++ b/source/blender/editors/mesh/editmesh_intersect.cc
@@ -652,7 +652,12 @@ static void bm_face_split_by_edges_island_connect(
       BMFace *f_pair[2];
       if (BM_edge_face_pair(edge_arr[i], &f_pair[0], &f_pair[1])) {
         if (BM_face_share_vert_count(f_pair[0], f_pair[1]) == 2) {
-          BMFace *f_new = BM_faces_join(bm, f_pair, 2, true);
+          BMFace *f_double;
+          BMFace *f_new = BM_faces_join(bm, f_pair, 2, true, &f_double);
+          /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+          BLI_assert_msg(f_double == nullptr,
+                         "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
           if (f_new) {
             BM_face_select_set(bm, f_new, true);
           }
diff --git a/source/blender/editors/mesh/editmesh_polybuild.cc b/source/blender/editors/mesh/editmesh_polybuild.cc
index e8898d591ef..288511dec7e 100644
--- a/source/blender/editors/mesh/editmesh_polybuild.cc
+++ b/source/blender/editors/mesh/editmesh_polybuild.cc
@@ -550,7 +550,12 @@ static wmOperatorStatus edbm_polybuild_dissolve_at_cursor_invoke(bContext *C,
     BMEdge *e_act = (BMEdge *)ele_act;
     BMLoop *l_a, *l_b;
     if (BM_edge_loop_pair(e_act, &l_a, &l_b)) {
-      BMFace *f_new = BM_faces_join_pair(bm, l_a, l_b, true);
+      BMFace *f_double;
+      BMFace *f_new = BM_faces_join_pair(bm, l_a, l_b, true, &f_double);
+      /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+      BLI_assert_msg(f_double == nullptr,
+                     "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
+
       if (f_new) {
         changed = true;
       }
diff --git a/source/blender/python/bmesh/bmesh_py_utils.cc b/source/blender/python/bmesh/bmesh_py_utils.cc
index 5434197ad29..f878ec7577b 100644
--- a/source/blender/python/bmesh/bmesh_py_utils.cc
+++ b/source/blender/python/bmesh/bmesh_py_utils.cc
@@ -653,7 +653,11 @@ static PyObject *bpy_bm_utils_face_join(PyObject * /*self*/, PyObject *args)
 
   /* Go ahead and join the face!
    * --------------------------- */
-  f_new = BM_faces_join(bm, face_array, int(face_seq_len), do_remove);
+  BMFace *f_double;
+  f_new = BM_faces_join(bm, face_array, int(face_seq_len), do_remove, &f_double);
+  /* See #BM_faces_join note on callers asserting when `r_double` is non-null. */
+  BLI_assert_msg(f_double == nullptr,
+                 "Doubled face detected at " AT ". Resulting mesh may be corrupt.");
 
   PyMem_FREE(face_array);