1 files changed, 288 insertions, 172 deletions

diff --git a/source/blender/bmesh/intern/bmesh_polygon.c b/source/blender/bmesh/intern/bmesh_polygon.c
index b0eddf73960..62b29e61d08 100644
--- a/source/blender/bmesh/intern/bmesh_polygon.c
+++ b/source/blender/bmesh/intern/bmesh_polygon.c
@@ -38,10 +38,13 @@
 #include "BLI_memarena.h"
 #include "BLI_polyfill2d.h"
 #include "BLI_polyfill2d_beautify.h"
+#include "BLI_linklist.h"
 
 #include "bmesh.h"
 #include "bmesh_tools.h"
 
+#include "BKE_customdata.h"
+
 #include "intern/bmesh_private.h"
 
 /**
@@ -103,7 +106,7 @@ static float bm_face_calc_poly_normal(const BMFace *f, float n[3])
  * but takes an array of vertex locations.
  */
 static float bm_face_calc_poly_normal_vertex_cos(
-        BMFace *f, float r_no[3],
+        const BMFace *f, float r_no[3],
         float const (*vertexCos)[3])
 {
 	BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
@@ -129,15 +132,15 @@ static float bm_face_calc_poly_normal_vertex_cos(
  * \brief COMPUTE POLY CENTER (BMFace)
  */
 static void bm_face_calc_poly_center_mean_vertex_cos(
-        BMFace *f, float r_cent[3],
+        const BMFace *f, float r_cent[3],
         float const (*vertexCos)[3])
 {
-	BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
-	BMLoop *l_iter  = l_first;
+	const BMLoop *l_first, *l_iter;
 
 	zero_v3(r_cent);
 
 	/* Newell's Method */
+	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
 	do {
 		add_v3_v3(r_cent, vertexCos[BM_elem_index_get(l_iter->v)]);
 	} while ((l_iter = l_iter->next) != l_first);
@@ -147,10 +150,14 @@ static void bm_face_calc_poly_center_mean_vertex_cos(
 /**
  * For tools that insist on using triangles, ideally we would cache this data.
  *
- * \param r_loops  Store face loop pointers, (f->len)
- * \param r_index  Store triangle triples, indices into \a r_loops,  ((f->len - 2) * 3)
+ * \param use_fixed_quad: When true, always split quad along (0 -> 2) regardless of concave corners,
+ * (as done in #BM_mesh_calc_tessellation).
+ * \param r_loops: Store face loop pointers, (f->len)
+ * \param r_index: Store triangle triples, indices into \a r_loops,  `((f->len - 2) * 3)`
  */
-void BM_face_calc_tessellation(const BMFace *f, BMLoop **r_loops, unsigned int (*r_index)[3])
+void BM_face_calc_tessellation(
+        const BMFace *f, const bool use_fixed_quad,
+        BMLoop **r_loops, unsigned int (*r_index)[3])
 {
 	BMLoop *l_first = BM_FACE_FIRST_LOOP(f);
 	BMLoop *l_iter;
@@ -164,7 +171,7 @@ void BM_face_calc_tessellation(const BMFace *f, BMLoop **r_loops, unsigned int (
 		r_index[0][1] = 1;
 		r_index[0][2] = 2;
 	}
-	else if (f->len == 4) {
+	else if (f->len == 4 && use_fixed_quad) {
 		*r_loops++ = (l_iter = l_first);
 		*r_loops++ = (l_iter = l_iter->next);
 		*r_loops++ = (l_iter = l_iter->next);
@@ -199,11 +206,51 @@ void BM_face_calc_tessellation(const BMFace *f, BMLoop **r_loops, unsigned int (
 }
 
 /**
+ * Return a point inside the face.
+ */
+void BM_face_calc_point_in_face(const BMFace *f, float r_co[3])
+{
+	const BMLoop *l_tri[3];
+
+	if (f->len == 3) {
+		const BMLoop *l = BM_FACE_FIRST_LOOP(f);
+		ARRAY_SET_ITEMS(l_tri, l, l->next, l->prev);
+	}
+	else {
+		/* tessellation here seems overkill when in many cases this will be the center,
+		 * but without this we can't be sure the point is inside a concave face. */
+		const int tottri = f->len - 2;
+		BMLoop **loops = BLI_array_alloca(loops, f->len);
+		unsigned int  (*index)[3] = BLI_array_alloca(index, tottri);
+		int j;
+		int j_best = 0;  /* use as fallback when unset */
+		float area_best  = -1.0f;
+
+		BM_face_calc_tessellation(f, false, loops, index);
+
+		for (j = 0; j < tottri; j++) {
+			const float *p1 = loops[index[j][0]]->v->co;
+			const float *p2 = loops[index[j][1]]->v->co;
+			const float *p3 = loops[index[j][2]]->v->co;
+			const float area = area_squared_tri_v3(p1, p2, p3);
+			if (area > area_best) {
+				j_best = j;
+				area_best = area;
+			}
+		}
+
+		ARRAY_SET_ITEMS(l_tri, loops[index[j_best][0]], loops[index[j_best][1]], loops[index[j_best][2]]);
+	}
+
+	mid_v3_v3v3v3(r_co, l_tri[0]->v->co, l_tri[1]->v->co, l_tri[2]->v->co);
+}
+
+/**
  * get the area of the face
  */
-float BM_face_calc_area(BMFace *f)
+float BM_face_calc_area(const BMFace *f)
 {
-	BMLoop *l_iter, *l_first;
+	const BMLoop *l_iter, *l_first;
 	float (*verts)[3] = BLI_array_alloca(verts, f->len);
 	float area;
 	unsigned int i = 0;
@@ -226,9 +273,9 @@ float BM_face_calc_area(BMFace *f)
 /**
  * compute the perimeter of an ngon
  */
-float BM_face_calc_perimeter(BMFace *f)
+float BM_face_calc_perimeter(const BMFace *f)
 {
-	BMLoop *l_iter, *l_first;
+	const BMLoop *l_iter, *l_first;
 	float perimeter = 0.0f;
 
 	l_iter = l_first = BM_FACE_FIRST_LOOP(f);
@@ -262,12 +309,12 @@ void BM_vert_tri_calc_plane(BMVert *verts[3], float r_plane[3])
  * Compute a meaningful direction along the face (use for manipulator axis).
  * \note result isnt normalized.
  */
-void BM_face_calc_plane(BMFace *f, float r_plane[3])
+void BM_face_calc_plane(const BMFace *f, float r_plane[3])
 {
 	if (f->len == 3) {
 		BMVert *verts[3];
 
-		BM_face_as_array_vert_tri(f, verts);
+		BM_face_as_array_vert_tri((BMFace *)f, verts);
 
 		BM_vert_tri_calc_plane(verts, r_plane);
 	}
@@ -276,7 +323,7 @@ void BM_face_calc_plane(BMFace *f, float r_plane[3])
 		float vec[3], vec_a[3], vec_b[3];
 
 		// BM_iter_as_array(NULL, BM_VERTS_OF_FACE, efa, (void **)verts, 4);
-		BM_face_as_array_vert_quad(f, verts);
+		BM_face_as_array_vert_quad((BMFace *)f, verts);
 
 		sub_v3_v3v3(vec_a, verts[3]->co, verts[2]->co);
 		sub_v3_v3v3(vec_b, verts[0]->co, verts[1]->co);
@@ -291,7 +338,7 @@ void BM_face_calc_plane(BMFace *f, float r_plane[3])
 		}
 	}
 	else {
-		BMLoop *l_long  = BM_face_find_longest_loop(f);
+		const BMLoop *l_long  = BM_face_find_longest_loop((BMFace *)f);
 
 		sub_v3_v3v3(r_plane, l_long->v->co, l_long->next->v->co);
 	}
@@ -302,10 +349,9 @@ void BM_face_calc_plane(BMFace *f, float r_plane[3])
 /**
  * computes center of face in 3d.  uses center of bounding box.
  */
-void BM_face_calc_center_bounds(BMFace *f, float r_cent[3])
+void BM_face_calc_center_bounds(const BMFace *f, float r_cent[3])
 {
-	BMLoop *l_iter;
-	BMLoop *l_first;
+	const BMLoop *l_iter, *l_first;
 	float min[3], max[3];
 
 	INIT_MINMAX(min, max);
@@ -321,9 +367,9 @@ void BM_face_calc_center_bounds(BMFace *f, float r_cent[3])
 /**
  * computes the center of a face, using the mean average
  */
-void BM_face_calc_center_mean(BMFace *f, float r_cent[3])
+void BM_face_calc_center_mean(const BMFace *f, float r_cent[3])
 {
-	BMLoop *l_iter, *l_first;
+	const BMLoop *l_iter, *l_first;
 
 	zero_v3(r_cent);
 
@@ -338,10 +384,10 @@ void BM_face_calc_center_mean(BMFace *f, float r_cent[3])
  * computes the center of a face, using the mean average
  * weighted by edge length
  */
-void BM_face_calc_center_mean_weighted(BMFace *f, float r_cent[3])
+void BM_face_calc_center_mean_weighted(const BMFace *f, float r_cent[3])
 {
-	BMLoop *l_iter;
-	BMLoop *l_first;
+	const BMLoop *l_iter;
+	const BMLoop *l_first;
 	float totw = 0.0f;
 	float w_prev;
 
@@ -422,32 +468,42 @@ void BM_edge_normals_update(BMEdge *e)
 	BM_vert_normal_update(e->v2);
 }
 
-bool BM_vert_normal_update_ex(BMVert *v, const char hflag, float r_no[3])
+static void bm_loop_normal_accum(const BMLoop *l, float no[3])
 {
-	/* TODO, we can normalize each edge only once, then compare with previous edge */
-
-	BMIter liter;
-	BMLoop *l;
-	int len = 0;
+	float vec1[3], vec2[3], fac;
 
-	zero_v3(r_no);
+	/* Same calculation used in BM_mesh_normals_update */
+	sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
+	sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
+	normalize_v3(vec1);
+	normalize_v3(vec2);
 
-	BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
-		if (BM_elem_flag_test(l->f, hflag)) {
-			float vec1[3], vec2[3], fac;
+	fac = saacos(-dot_v3v3(vec1, vec2));
 
-			/* Same calculation used in BM_mesh_normals_update */
-			sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
-			sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
-			normalize_v3(vec1);
-			normalize_v3(vec2);
+	madd_v3_v3fl(no, l->f->no, fac);
+}
 
-			fac = saacos(-dot_v3v3(vec1, vec2));
+bool BM_vert_calc_normal_ex(const BMVert *v, const char hflag, float r_no[3])
+{
+	int len = 0;
 
-			madd_v3_v3fl(r_no, l->f->no, fac);
+	zero_v3(r_no);
 
-			len++;
-		}
+	if (v->e) {
+		const BMEdge *e = v->e;
+		do {
+			if (e->l) {
+				const BMLoop *l = e->l;
+				do {
+					if (l->v == v) {
+						if (BM_elem_flag_test(l->f, hflag)) {
+							bm_loop_normal_accum(l, r_no);
+							len++;
+						}
+					}
+				} while ((l = l->radial_next) != e->l);
+			}
+		} while ((e = bmesh_disk_edge_next(e, v)) != v->e);
 	}
 
 	if (len) {
@@ -459,33 +515,56 @@ bool BM_vert_normal_update_ex(BMVert *v, const char hflag, float r_no[3])
 	}
 }
 
-/**
- * update a vert normal (but not the faces incident on it)
- */
-void BM_vert_normal_update(BMVert *v)
+bool BM_vert_calc_normal(const BMVert *v, float r_no[3])
 {
-	/* TODO, we can normalize each edge only once, then compare with previous edge */
-
-	BMIter liter;
-	BMLoop *l;
 	int len = 0;
 
-	zero_v3(v->no);
+	zero_v3(r_no);
 
-	BM_ITER_ELEM (l, &liter, v, BM_LOOPS_OF_VERT) {
-		float vec1[3], vec2[3], fac;
+	if (v->e) {
+		const BMEdge *e = v->e;
+		do {
+			if (e->l) {
+				const BMLoop *l = e->l;
+				do {
+					if (l->v == v) {
+						bm_loop_normal_accum(l, r_no);
+						len++;
+					}
+				} while ((l = l->radial_next) != e->l);
+			}
+		} while ((e = bmesh_disk_edge_next(e, v)) != v->e);
+	}
 
-		/* Same calculation used in BM_mesh_normals_update */
-		sub_v3_v3v3(vec1, l->v->co, l->prev->v->co);
-		sub_v3_v3v3(vec2, l->next->v->co, l->v->co);
-		normalize_v3(vec1);
-		normalize_v3(vec2);
+	if (len) {
+		normalize_v3(r_no);
+		return true;
+	}
+	else {
+		return false;
+	}
+}
 
-		fac = saacos(-dot_v3v3(vec1, vec2));
+void BM_vert_normal_update_all(BMVert *v)
+{
+	int len = 0;
 
-		madd_v3_v3fl(v->no, l->f->no, fac);
+	zero_v3(v->no);
 
-		len++;
+	if (v->e) {
+		const BMEdge *e = v->e;
+		do {
+			if (e->l) {
+				const BMLoop *l = e->l;
+				do {
+					if (l->v == v) {
+						BM_face_normal_update(l->f);
+						bm_loop_normal_accum(l, v->no);
+						len++;
+					}
+				} while ((l = l->radial_next) != e->l);
+			}
+		} while ((e = bmesh_disk_edge_next(e, v)) != v->e);
 	}
 
 	if (len) {
@@ -493,16 +572,12 @@ void BM_vert_normal_update(BMVert *v)
 	}
 }
 
-void BM_vert_normal_update_all(BMVert *v)
+/**
+ * update a vert normal (but not the faces incident on it)
+ */
+void BM_vert_normal_update(BMVert *v)
 {
-	BMIter iter;
-	BMFace *f;
-
-	BM_ITER_ELEM (f, &iter, v, BM_FACES_OF_VERT) {
-		BM_face_normal_update(f);
-	}
-
-	BM_vert_normal_update(v);
+	BM_vert_calc_normal(v, v->no);
 }
 
 /**
@@ -551,7 +626,7 @@ void BM_face_normal_update(BMFace *f)
 
 /* exact same as 'BM_face_calc_normal' but accepts vertex coords */
 float BM_face_calc_normal_vcos(
-        BMesh *bm, BMFace *f, float r_no[3],
+        const BMesh *bm, const BMFace *f, float r_no[3],
         float const (*vertexCos)[3])
 {
 	BMLoop *l;
@@ -589,13 +664,13 @@ float BM_face_calc_normal_vcos(
 /**
  * Calculates the face subset normal.
  */
-float BM_face_calc_normal_subset(BMLoop *l_first, BMLoop *l_last, float r_no[3])
+float BM_face_calc_normal_subset(const BMLoop *l_first, const BMLoop *l_last, float r_no[3])
 {
 	const float *v_prev, *v_curr;
 
 	/* Newell's Method */
-	BMLoop *l_iter = l_first;
-	BMLoop *l_term = l_last->next;
+	const BMLoop *l_iter = l_first;
+	const BMLoop *l_term = l_last->next;
 
 	zero_v3(r_no);
 
@@ -611,7 +686,7 @@ float BM_face_calc_normal_subset(BMLoop *l_first, BMLoop *l_last, float r_no[3])
 
 /* exact same as 'BM_face_calc_normal' but accepts vertex coords */
 void BM_face_calc_center_mean_vcos(
-        BMesh *bm, BMFace *f, float r_cent[3],
+        const BMesh *bm, const BMFace *f, float r_cent[3],
         float const (*vertexCos)[3])
 {
 	/* must have valid index data */
@@ -627,12 +702,20 @@ void BM_face_calc_center_mean_vcos(
  * Reverses the winding of a face.
  * \note This updates the calculated normal.
  */
-void BM_face_normal_flip(BMesh *bm, BMFace *f)
+void BM_face_normal_flip_ex(
+        BMesh *bm, BMFace *f,
+        const int cd_loop_mdisp_offset, const bool use_loop_mdisp_flip)
 {
-	bmesh_loop_reverse(bm, f);
+	bmesh_loop_reverse(bm, f, cd_loop_mdisp_offset, use_loop_mdisp_flip);
 	negate_v3(f->no);
 }
 
+void BM_face_normal_flip(BMesh *bm, BMFace *f)
+{
+	const int cd_loop_mdisp_offset = CustomData_get_offset(&bm->ldata, CD_MDISPS);
+	BM_face_normal_flip_ex(bm, f, cd_loop_mdisp_offset, true);
+}
+
 /* detects if two line segments cross each other (intersects).
  * note, there could be more winding cases then there needs to be. */
 static bool line_crosses_v2f(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
@@ -704,7 +787,7 @@ static bool line_crosses_v2f(const float v1[2], const float v2[2], const float v
  * instead of projecting co directly into f's orientation space,
  * so there might be accuracy issues.
  */
-bool BM_face_point_inside_test(BMFace *f, const float co[3])
+bool BM_face_point_inside_test(const BMFace *f, const float co[3])
 {
 	float axis_mat[3][3];
 	float (*projverts)[2] = BLI_array_alloca(projverts, f->len);
@@ -713,8 +796,7 @@ bool BM_face_point_inside_test(BMFace *f, const float co[3])
 	BMLoop *l_iter;
 	int i;
 	
-	if (is_zero_v3(f->no))
-		BM_face_normal_update(f);
+	BLI_assert(BM_face_is_normal_valid(f));
 
 	axis_dominant_v3_to_m3(axis_mat, f->no);
 
@@ -738,6 +820,12 @@ bool BM_face_point_inside_test(BMFace *f, const float co[3])
  * with a length equal to (f->len - 3). It will be filled with the new
  * triangles (not including the original triangle).
  *
+ * \param r_faces_double: When newly created faces are duplicates of existing faces, they're added to this list.
+ * Caller must handle de-duplication.
+ * This is done because its possible _all_ faces exist already,
+ * and in that case we would have to remove all faces including the one passed,
+ * which causes complications adding/removing faces while looking over them.
+ *
  * \note The number of faces is _almost_ always (f->len - 3),
  *       However there may be faces that already occupying the
  *       triangles we would make, so the caller must check \a r_faces_new_tot.
@@ -750,6 +838,7 @@ void BM_face_triangulate(
         int     *r_faces_new_tot,
         BMEdge **r_edges_new,
         int     *r_edges_new_tot,
+        LinkNode **r_faces_double,
         const int quad_method,
         const int ngon_method,
         const bool use_tag,
@@ -758,125 +847,132 @@ void BM_face_triangulate(
         /* use for MOD_TRIANGULATE_NGON_BEAUTY only! */
         struct Heap *pf_heap, struct EdgeHash *pf_ehash)
 {
-	BMLoop *l_iter, *l_first, *l_new;
+	const int cd_loop_mdisp_offset = CustomData_get_offset(&bm->ldata, CD_MDISPS);
+	const bool use_beauty = (ngon_method == MOD_TRIANGULATE_NGON_BEAUTY);
+	BMLoop *l_first, *l_new;
 	BMFace *f_new;
 	int nf_i = 0;
 	int ne_i = 0;
-	bool use_beauty = (ngon_method == MOD_TRIANGULATE_NGON_BEAUTY);
 
 	BLI_assert(BM_face_is_normal_valid(f));
 
 	/* ensure both are valid or NULL */
 	BLI_assert((r_faces_new == NULL) == (r_faces_new_tot == NULL));
 
-	if (f->len == 4) {
-		BMLoop *l_v1, *l_v2;
-		l_first = BM_FACE_FIRST_LOOP(f);
+	BLI_assert(f->len > 3);
 
-		switch (quad_method) {
-			case MOD_TRIANGULATE_QUAD_FIXED:
-			{
-				l_v1 = l_first;
-				l_v2 = l_first->next->next;
-				break;
-			}
-			case MOD_TRIANGULATE_QUAD_ALTERNATE:
-			{
-				l_v1 = l_first->next;
-				l_v2 = l_first->prev;
-				break;
-			}
-			case MOD_TRIANGULATE_QUAD_SHORTEDGE:
-			case MOD_TRIANGULATE_QUAD_BEAUTY:
-			default:
-			{
-				BMLoop *l_v3, *l_v4;
-				bool split_24;
-
-				l_v1 = l_first->next;
-				l_v2 = l_first->next->next;
-				l_v3 = l_first->prev;
-				l_v4 = l_first;
-
-				if (quad_method == MOD_TRIANGULATE_QUAD_SHORTEDGE) {
-					float d1, d2;
-					d1 = len_squared_v3v3(l_v4->v->co, l_v2->v->co);
-					d2 = len_squared_v3v3(l_v1->v->co, l_v3->v->co);
-					split_24 = ((d2 - d1) > 0.0f);
+	{
+		BMLoop **loops = BLI_array_alloca(loops, f->len);
+		unsigned int (*tris)[3] = BLI_array_alloca(tris, f->len);
+		const int totfilltri = f->len - 2;
+		const int last_tri = f->len - 3;
+		int i;
+		/* for mdisps */
+		float f_center[3];
+
+		if (f->len == 4) {
+			/* even though we're not using BLI_polyfill, fill in 'tris' and 'loops'
+			 * so we can share code to handle face creation afterwards. */
+			BMLoop *l_v1, *l_v2;
+
+			l_first = BM_FACE_FIRST_LOOP(f);
+
+			switch (quad_method) {
+				case MOD_TRIANGULATE_QUAD_FIXED:
+				{
+					l_v1 = l_first;
+					l_v2 = l_first->next->next;
+					break;
 				}
-				else {
-					/* first check if the quad is concave on either diagonal */
-					const int flip_flag = is_quad_flip_v3(l_v1->v->co, l_v2->v->co, l_v3->v->co, l_v4->v->co);
-					if (UNLIKELY(flip_flag & (1 << 0))) {
-						split_24 = true;
-					}
-					else if (UNLIKELY(flip_flag & (1 << 1))) {
-						split_24 = false;
+				case MOD_TRIANGULATE_QUAD_ALTERNATE:
+				{
+					l_v1 = l_first->next;
+					l_v2 = l_first->prev;
+					break;
+				}
+				case MOD_TRIANGULATE_QUAD_SHORTEDGE:
+				case MOD_TRIANGULATE_QUAD_BEAUTY:
+				default:
+				{
+					BMLoop *l_v3, *l_v4;
+					bool split_24;
+
+					l_v1 = l_first->next;
+					l_v2 = l_first->next->next;
+					l_v3 = l_first->prev;
+					l_v4 = l_first;
+
+					if (quad_method == MOD_TRIANGULATE_QUAD_SHORTEDGE) {
+						float d1, d2;
+						d1 = len_squared_v3v3(l_v4->v->co, l_v2->v->co);
+						d2 = len_squared_v3v3(l_v1->v->co, l_v3->v->co);
+						split_24 = ((d2 - d1) > 0.0f);
 					}
 					else {
-						split_24 = (BM_verts_calc_rotate_beauty(l_v1->v, l_v2->v, l_v3->v, l_v4->v, 0, 0) > 0.0f);
+						/* first check if the quad is concave on either diagonal */
+						const int flip_flag = is_quad_flip_v3(l_v1->v->co, l_v2->v->co, l_v3->v->co, l_v4->v->co);
+						if (UNLIKELY(flip_flag & (1 << 0))) {
+							split_24 = true;
+						}
+						else if (UNLIKELY(flip_flag & (1 << 1))) {
+							split_24 = false;
+						}
+						else {
+							split_24 = (BM_verts_calc_rotate_beauty(l_v1->v, l_v2->v, l_v3->v, l_v4->v, 0, 0) > 0.0f);
+						}
 					}
-				}
 
-				/* named confusingly, l_v1 is in fact the second vertex */
-				if (split_24) {
-					l_v1 = l_v4;
-					//l_v2 = l_v2;
-				}
-				else {
-					//l_v1 = l_v1;
-					l_v2 = l_v3;
+					/* named confusingly, l_v1 is in fact the second vertex */
+					if (split_24) {
+						l_v1 = l_v4;
+						//l_v2 = l_v2;
+					}
+					else {
+						//l_v1 = l_v1;
+						l_v2 = l_v3;
+					}
+					break;
 				}
-				break;
 			}
-		}
-
-		f_new = BM_face_split(bm, f, l_v1, l_v2, &l_new, NULL, true);
-		copy_v3_v3(f_new->no, f->no);
 
-		if (use_tag) {
-			BM_elem_flag_enable(l_new->e, BM_ELEM_TAG);
-			BM_elem_flag_enable(f_new, BM_ELEM_TAG);
-		}
+			loops[0] = l_v1;
+			loops[1] = l_v1->next;
+			loops[2] = l_v2;
+			loops[3] = l_v2->next;
 
-		if (r_faces_new) {
-			r_faces_new[nf_i++] = f_new;
-		}
-		if (r_edges_new) {
-			r_edges_new[ne_i++] = l_new->e;
+			ARRAY_SET_ITEMS(tris[0], 0, 1, 2);
+			ARRAY_SET_ITEMS(tris[1], 0, 2, 3);
 		}
-	}
-	else if (f->len > 4) {
+		else {
+			BMLoop *l_iter;
+			float axis_mat[3][3];
+			float (*projverts)[2] = BLI_array_alloca(projverts, f->len);
 
-		float axis_mat[3][3];
-		float (*projverts)[2] = BLI_array_alloca(projverts, f->len);
-		BMLoop **loops = BLI_array_alloca(loops, f->len);
-		unsigned int (*tris)[3] = BLI_array_alloca(tris, f->len);
-		const int totfilltri = f->len - 2;
-		const int last_tri = f->len - 3;
-		int i;
+			axis_dominant_v3_to_m3_negate(axis_mat, f->no);
 
-		axis_dominant_v3_to_m3_negate(axis_mat, f->no);
+			for (i = 0, l_iter = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l_iter = l_iter->next) {
+				loops[i] = l_iter;
+				mul_v2_m3v3(projverts[i], axis_mat, l_iter->v->co);
+			}
 
-		for (i = 0, l_iter = BM_FACE_FIRST_LOOP(f); i < f->len; i++, l_iter = l_iter->next) {
-			loops[i] = l_iter;
-			mul_v2_m3v3(projverts[i], axis_mat, l_iter->v->co);
-		}
+			BLI_polyfill_calc_arena((const float (*)[2])projverts, f->len, 1, tris,
+			                        pf_arena);
 
-		BLI_polyfill_calc_arena((const float (*)[2])projverts, f->len, 1, tris,
-		                        pf_arena);
+			if (use_beauty) {
+				BLI_polyfill_beautify(
+				        (const float (*)[2])projverts, f->len, tris,
+				        pf_arena, pf_heap, pf_ehash);
+			}
 
-		if (use_beauty) {
-			BLI_polyfill_beautify(
-			        (const float (*)[2])projverts, f->len, tris,
-			        pf_arena, pf_heap, pf_ehash);
+			BLI_memarena_clear(pf_arena);
 		}
 
-		BLI_memarena_clear(pf_arena);
+		if (cd_loop_mdisp_offset != -1) {
+			BM_face_calc_center_mean(f, f_center);
+		}
 
 		/* loop over calculated triangles and create new geometry */
 		for (i = 0; i < totfilltri; i++) {
-			/* the order is reverse, otherwise the normal is flipped */
 			BMLoop *l_tri[3] = {
 			    loops[tris[i][0]],
 			    loops[tris[i][1]],
@@ -892,6 +988,18 @@ void BM_face_triangulate(
 
 			BLI_assert(v_tri[0] == l_new->v);
 
+			/* check for duplicate */
+			if (l_new->radial_next != l_new) {
+				BMLoop *l_iter = l_new->radial_next;
+				do {
+					if (UNLIKELY((l_iter->f->len == 3) && (l_new->prev->v == l_iter->prev->v))) {
+						/* Check the last tri because we swap last f_new with f at the end... */
+						BLI_linklist_prepend(r_faces_double, (i != last_tri) ? f_new : f);
+						break;
+					}
+				} while ((l_iter = l_iter->radial_next) != l_new);
+			}
+
 			/* copy CD data */
 			BM_elem_attrs_copy(bm, bm, l_tri[0], l_new);
 			BM_elem_attrs_copy(bm, bm, l_tri[1], l_new->next);
@@ -909,6 +1017,8 @@ void BM_face_triangulate(
 
 			if (use_tag || r_edges_new) {
 				/* new faces loops */
+				BMLoop *l_iter;
+
 				l_iter = l_first = l_new;
 				do {
 					BMEdge *e = l_iter->e;
@@ -927,6 +1037,12 @@ void BM_face_triangulate(
 					/* note, never disable tag's */
 				} while ((l_iter = l_iter->next) != l_first);
 			}
+
+			if (cd_loop_mdisp_offset != -1) {
+				float f_new_center[3];
+				BM_face_calc_center_mean(f_new, f_new_center);
+				BM_face_interp_multires_ex(bm, f_new, f, f_new_center, f_center, cd_loop_mdisp_offset);
+			}
 		}
 
 		{
@@ -1159,12 +1275,12 @@ void BM_face_as_array_loop_quad(BMFace *f, BMLoop *r_loops[4])
 
 
 /**
- * \brief BM_bmesh_calc_tessellation get the looptris and its number from a certain bmesh
+ * \brief BM_mesh_calc_tessellation get the looptris and its number from a certain bmesh
  * \param looptris
  *
  * \note \a looptris  Must be pre-allocated to at least the size of given by: poly_to_tri_count
  */
-void BM_bmesh_calc_tessellation(BMesh *bm, BMLoop *(*looptris)[3], int *r_looptris_tot)
+void BM_mesh_calc_tessellation(BMesh *bm, BMLoop *(*looptris)[3], int *r_looptris_tot)
 {
 	/* use this to avoid locking pthread for _every_ polygon
 	 * and calling the fill function */