Add Custom Loop Normals.

This is the core code for it, tools (datatransfer and modifier) will come in next commits).
RNA api is already there, though.

See the code for details, but basically, we define, for each 'smooth fan'
(which is a set of adjacent loops around a same vertex that are smooth, i.e. have a single same normal),
a 'loop normal space' (or lnor space), using auto-computed normal and relevant edges, and store
custom normal as two angular factors inside that space. This allows to have custom normals
'following' deformations of the geometry, and to only save two shorts per loop in new clnor CDLayer.

Normal manipulation (editing, mixing, interpolating, etc.) shall always happen with plain 3D vectors normals,
and be converted back into storage format at the end.

Clnor computation has also been threaded (at least for Mesh case, not for BMesh), since the process can
be rather heavy with high poly meshes.

Also, bumping subversion, and fix mess in 2.70 versioning code.

21 files changed, 1684 insertions, 304 deletions

diff --git a/source/blender/blenkernel/BKE_DerivedMesh.h b/source/blender/blenkernel/BKE_DerivedMesh.h
index da982f869e8..04251641d9b 100644
--- a/source/blender/blenkernel/BKE_DerivedMesh.h
+++ b/source/blender/blenkernel/BKE_DerivedMesh.h
@@ -88,6 +88,7 @@ struct MTFace;
 struct Object;
 struct Scene;
 struct Mesh;
+struct MLoopNorSpaceArray;
 struct BMEditMesh;
 struct KeyBlock;
 struct ModifierData;
@@ -199,6 +200,10 @@ struct DerivedMesh {
 	/** Calculate loop (split) normals */
 	void (*calcLoopNormals)(DerivedMesh *dm, const bool use_split_normals, const float split_angle);
 
+	/** Calculate loop (split) normals, and returns split loop normal spacearr. */
+	void (*calcLoopNormalsSpaceArray)(DerivedMesh *dm, const bool use_split_normals, const float split_angle,
+	                                  struct MLoopNorSpaceArray *r_lnors_spacearr);
+
 	/** Recalculates mesh tessellation */
 	void (*recalcTessellation)(DerivedMesh *dm);
 
diff --git a/source/blender/blenkernel/BKE_blender.h b/source/blender/blenkernel/BKE_blender.h
index a2214f79618..0c10b209fb7 100644
--- a/source/blender/blenkernel/BKE_blender.h
+++ b/source/blender/blenkernel/BKE_blender.h
@@ -42,7 +42,7 @@ extern "C" {
  * and keep comment above the defines.
  * Use STRINGIFY() rather than defining with quotes */
 #define BLENDER_VERSION         273
-#define BLENDER_SUBVERSION      5
+#define BLENDER_SUBVERSION      6
 /* 262 was the last editmesh release but it has compatibility code for bmesh data */
 #define BLENDER_MINVERSION      270
 #define BLENDER_MINSUBVERSION   5
diff --git a/source/blender/blenkernel/BKE_cdderivedmesh.h b/source/blender/blenkernel/BKE_cdderivedmesh.h
index 448617f73f9..59ec316a403 100644
--- a/source/blender/blenkernel/BKE_cdderivedmesh.h
+++ b/source/blender/blenkernel/BKE_cdderivedmesh.h
@@ -40,6 +40,7 @@
 struct DerivedMesh;
 struct BMEditMesh;
 struct Mesh;
+struct MLoopNorSpaceArray;
 struct Object;
 
 /* creates a new CDDerivedMesh */
@@ -107,6 +108,8 @@ void CDDM_calc_normals(struct DerivedMesh *dm);
 void CDDM_calc_normals_tessface(struct DerivedMesh *dm);
 
 void CDDM_calc_loop_normals(struct DerivedMesh *dm, const bool use_split_normals, const float split_angle);
+void CDDM_calc_loop_normals_spacearr(struct DerivedMesh *dm, const bool use_split_normals, const float split_angle,
+                                     struct MLoopNorSpaceArray *r_lnors_spacearr);
 
 /* calculates edges for a CDDerivedMesh (from face data)
  * this completely replaces the current edge data in the DerivedMesh
diff --git a/source/blender/blenkernel/BKE_mesh.h b/source/blender/blenkernel/BKE_mesh.h
index 3138fd69cc4..224be0f3685 100644
--- a/source/blender/blenkernel/BKE_mesh.h
+++ b/source/blender/blenkernel/BKE_mesh.h
@@ -36,6 +36,9 @@ struct BoundBox;
 struct DispList;
 struct EdgeHash;
 struct ListBase;
+struct LinkNode;
+struct BLI_Stack;
+struct MemArena;
 struct BMEditMesh;
 struct BMesh;
 struct Main;
@@ -174,11 +177,6 @@ void BKE_mesh_calc_normals_tessface(
         struct MVert *mverts, int numVerts,
         struct MFace *mfaces, int numFaces,
         float (*r_faceNors)[3]);
-void BKE_mesh_normals_loop_split(
-        struct MVert *mverts, const int numVerts, struct MEdge *medges, const int numEdges,
-        struct MLoop *mloops, float (*r_loopnors)[3], const int numLoops,
-        struct MPoly *mpolys, float (*polynors)[3], const int numPolys,
-        const bool use_split_normals, float split_angle);
 void BKE_mesh_loop_tangents_ex(
         struct MVert *mverts, const int numVerts, struct MLoop *mloops, float (*r_looptangent)[4], float (*loopnors)[3],
         struct MLoopUV *loopuv, const int numLoops, struct MPoly *mpolys, const int numPolys,
@@ -186,6 +184,56 @@ void BKE_mesh_loop_tangents_ex(
 void BKE_mesh_loop_tangents(
         struct Mesh *mesh, const char *uvmap, float (*r_looptangents)[4], struct ReportList *reports);
 
+/**
+ * References a contiguous loop-fan with normal offset vars.
+ */
+typedef struct MLoopNorSpace {
+	float vec_lnor[3];      /* Automatically computed loop normal. */
+	float vec_ref[3];       /* Reference vector, orthogonal to vec_lnor. */
+	float vec_ortho[3];     /* Third vector, orthogonal to vec_lnor and vec_ref. */
+	float ref_alpha;        /* Reference angle, around vec_ortho, in ]0, pi] range (0.0 marks that space as invalid). */
+	float ref_beta;         /* Reference angle, around vec_lnor, in ]0, 2pi] range (0.0 marks that space as invalid). */
+	struct LinkNode *loops; /* All indices (uint_in_ptr) of loops using this lnor space (i.e. smooth fan of loops). */
+} MLoopNorSpace;
+/**
+ * Collection of #MLoopNorSpace basic storage & pre-allocation.
+ */
+typedef struct MLoopNorSpaceArray {
+	MLoopNorSpace **lspacearr;    /* MLoop aligned array */
+	struct LinkNode *loops_pool;  /* Allocated once, avoids to call BLI_linklist_prepend_arena() for each loop! */
+	struct MemArena *mem;
+} MLoopNorSpaceArray;
+void BKE_lnor_spacearr_init(MLoopNorSpaceArray *lnors_spacearr, const int numLoops);
+void BKE_lnor_spacearr_clear(MLoopNorSpaceArray *lnors_spacearr);
+void BKE_lnor_spacearr_free(MLoopNorSpaceArray *lnors_spacearr);
+MLoopNorSpace *BKE_lnor_space_create(MLoopNorSpaceArray *lnors_spacearr);
+void BKE_lnor_space_define(
+        MLoopNorSpace *lnor_space, const float lnor[3], float vec_ref[3], float vec_other[3],
+        struct BLI_Stack *edge_vectors);
+void BKE_lnor_space_add_loop(
+        MLoopNorSpaceArray *lnors_spacearr, MLoopNorSpace *lnor_space, const int ml_index, const bool add_to_list);
+void BKE_lnor_space_custom_data_to_normal(MLoopNorSpace *lnor_space, const short clnor_data[2], float r_custom_lnor[3]);
+void BKE_lnor_space_custom_normal_to_data(MLoopNorSpace *lnor_space, const float custom_lnor[3], short r_clnor_data[2]);
+
+bool BKE_mesh_has_custom_loop_normals(struct Mesh *me);
+
+void BKE_mesh_normals_loop_split(struct MVert *mverts, const int numVerts, struct MEdge *medges, const int numEdges,
+        struct MLoop *mloops, float (*r_loopnors)[3], const int numLoops,
+        struct MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        const bool use_split_normals, float split_angle,
+        MLoopNorSpaceArray *r_lnors_spacearr, short (*clnors_data)[2], int *r_loop_to_poly);
+
+void BKE_mesh_normals_loop_custom_set(
+        struct MVert *mverts, const int numVerts, struct MEdge *medges, const int numEdges,
+        struct MLoop *mloops, float (*custom_loopnors)[3], const int numLoops,
+        struct MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        short (*r_clnors_data)[2]);
+void BKE_mesh_normals_loop_custom_from_vertices_set(
+        struct MVert *mverts, float (*custom_vertnors)[3], const int numVerts,
+        struct MEdge *medges, const int numEdges, struct MLoop *mloops, const int numLoops,
+        struct MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        short (*r_clnors_data)[2]);
+
 void BKE_mesh_calc_poly_normal(
         struct MPoly *mpoly, struct MLoop *loopstart,
         struct MVert *mvarray, float no[3]);
diff --git a/source/blender/blenkernel/intern/cdderivedmesh.c b/source/blender/blenkernel/intern/cdderivedmesh.c
index 084376e768e..2554151c99d 100644
--- a/source/blender/blenkernel/intern/cdderivedmesh.c
+++ b/source/blender/blenkernel/intern/cdderivedmesh.c
@@ -1540,6 +1540,7 @@ static CDDerivedMesh *cdDM_create(const char *desc)
 
 	dm->calcNormals = CDDM_calc_normals;
 	dm->calcLoopNormals = CDDM_calc_loop_normals;
+	dm->calcLoopNormalsSpaceArray = CDDM_calc_loop_normals_spacearr;
 	dm->recalcTessellation = CDDM_recalc_tessellation;
 
 	dm->getVertCos = cdDM_getVertCos;
@@ -2158,6 +2159,14 @@ void CDDM_calc_normals(DerivedMesh *dm)
 
 void CDDM_calc_loop_normals(DerivedMesh *dm, const bool use_split_normals, const float split_angle)
 {
+	CDDM_calc_loop_normals_spacearr(dm, use_split_normals, split_angle, NULL);
+}
+
+/* #define DEBUG_CLNORS */
+
+void CDDM_calc_loop_normals_spacearr(
+        DerivedMesh *dm, const bool use_split_normals, const float split_angle, MLoopNorSpaceArray *r_lnors_spacearr)
+{
 	MVert *mverts = dm->getVertArray(dm);
 	MEdge *medges = dm->getEdgeArray(dm);
 	MLoop *mloops = dm->getLoopArray(dm);
@@ -2166,6 +2175,7 @@ void CDDM_calc_loop_normals(DerivedMesh *dm, const bool use_split_normals, const
 	CustomData *ldata, *pdata;
 
 	float (*lnors)[3];
+	short (*clnor_data)[2];
 	float (*pnors)[3];
 
 	const int numVerts = dm->getNumVerts(dm);
@@ -2193,8 +2203,37 @@ void CDDM_calc_loop_normals(DerivedMesh *dm, const bool use_split_normals, const
 
 	dm->dirty &= ~DM_DIRTY_NORMALS;
 
+	clnor_data = CustomData_get_layer(ldata, CD_CUSTOMLOOPNORMAL);
+
 	BKE_mesh_normals_loop_split(mverts, numVerts, medges, numEdges, mloops, lnors, numLoops,
-	                            mpolys, pnors, numPolys, use_split_normals, split_angle);
+	                            mpolys, (const float (*)[3])pnors, numPolys,
+	                            use_split_normals, split_angle,
+	                            r_lnors_spacearr, clnor_data, NULL);
+#ifdef DEBUG_CLNORS
+	if (r_lnors_spacearr) {
+		int i;
+		for (i = 0; i < numLoops; i++) {
+			if (r_lnors_spacearr->lspacearr[i]->ref_alpha != 0.0f) {
+				LinkNode *loops = r_lnors_spacearr->lspacearr[i]->loops;
+				printf("Loop %d uses lnor space %p:\n", i, r_lnors_spacearr->lspacearr[i]);
+				print_v3("\tfinal lnor", lnors[i]);
+				print_v3("\tauto lnor", r_lnors_spacearr->lspacearr[i]->vec_lnor);
+				print_v3("\tref_vec", r_lnors_spacearr->lspacearr[i]->vec_ref);
+				printf("\talpha: %f\n\tbeta: %f\n\tloops: %p\n", r_lnors_spacearr->lspacearr[i]->ref_alpha,
+				       r_lnors_spacearr->lspacearr[i]->ref_beta, r_lnors_spacearr->lspacearr[i]->loops);
+				printf("\t\t(shared with loops");
+				while (loops) {
+					printf(" %d", GET_INT_FROM_POINTER(loops->link));
+					loops = loops->next;
+				}
+				printf(")\n");
+			}
+			else {
+				printf("Loop %d has no lnor space\n", i);
+			}
+		}
+	}
+#endif
 }
 
 
diff --git a/source/blender/blenkernel/intern/customdata.c b/source/blender/blenkernel/intern/customdata.c
index a259958dc0c..0c437f6e02f 100644
--- a/source/blender/blenkernel/intern/customdata.c
+++ b/source/blender/blenkernel/intern/customdata.c
@@ -1255,6 +1255,8 @@ static const LayerTypeInfo LAYERTYPEINFO[CD_NUMTYPES] = {
 	{sizeof(float[4]), "", 0, NULL, NULL, NULL, NULL, NULL, NULL},
 	/* 40: CD_TESSLOOPNORMAL */
 	{sizeof(short[4][3]), "", 0, NULL, NULL, NULL, NULL, layerSwap_flnor, NULL},
+	/* 41: CD_CUSTOMLOOPNORMAL */
+	{sizeof(short[2]), "vec2s", 1, NULL, NULL, NULL, NULL, NULL, NULL},
 };
 
 /* note, numbers are from trunk and need updating for bmesh */
@@ -1270,7 +1272,8 @@ static const char *LAYERTYPENAMES[CD_NUMTYPES] = {
 	/* 25-29 */ "CDMPoly", "CDMLoop", "CDShapeKeyIndex", "CDShapeKey", "CDBevelWeight",
 	/* 30-34 */ "CDSubSurfCrease", "CDOrigSpaceLoop", "CDPreviewLoopCol", "CDBMElemPyPtr", "CDPaintMask",
 	/* 35-36 */ "CDGridPaintMask", "CDMVertSkin",
-	/* 37-40 */ "CDFreestyleEdge", "CDFreestyleFace", "CDMLoopTangent", "CDTessLoopNormal",
+	/* 37-38 */ "CDFreestyleEdge", "CDFreestyleFace",
+	/* 39-41 */ "CDMLoopTangent", "CDTessLoopNormal", "CDCustomLoopNormal",
 };
 
 
@@ -1282,26 +1285,29 @@ const CustomDataMask CD_MASK_MESH =
     CD_MASK_PROP_FLT | CD_MASK_PROP_INT | CD_MASK_PROP_STR | CD_MASK_MDISPS |
     CD_MASK_MLOOPUV | CD_MASK_MLOOPCOL | CD_MASK_MPOLY | CD_MASK_MLOOP |
     CD_MASK_MTEXPOLY | CD_MASK_RECAST | CD_MASK_PAINT_MASK |
-    CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN | CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE;
+    CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN | CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE |
+    CD_MASK_CUSTOMLOOPNORMAL;
 const CustomDataMask CD_MASK_EDITMESH =
     CD_MASK_MSTICKY | CD_MASK_MDEFORMVERT | CD_MASK_MTFACE | CD_MASK_MLOOPUV |
     CD_MASK_MLOOPCOL | CD_MASK_MTEXPOLY | CD_MASK_SHAPE_KEYINDEX |
     CD_MASK_MCOL | CD_MASK_PROP_FLT | CD_MASK_PROP_INT | CD_MASK_PROP_STR |
     CD_MASK_MDISPS | CD_MASK_SHAPEKEY | CD_MASK_RECAST | CD_MASK_PAINT_MASK |
-    CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN;
+    CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN | CD_MASK_CUSTOMLOOPNORMAL;
 const CustomDataMask CD_MASK_DERIVEDMESH =
     CD_MASK_MSTICKY | CD_MASK_MDEFORMVERT | CD_MASK_MTFACE |
     CD_MASK_MCOL | CD_MASK_PROP_FLT | CD_MASK_PROP_INT | CD_MASK_CLOTH_ORCO |
     CD_MASK_MLOOPUV | CD_MASK_MLOOPCOL | CD_MASK_MTEXPOLY | CD_MASK_PREVIEW_MLOOPCOL |
     CD_MASK_PROP_STR | CD_MASK_ORIGSPACE | CD_MASK_ORIGSPACE_MLOOP | CD_MASK_ORCO | CD_MASK_TANGENT |
     CD_MASK_PREVIEW_MCOL | CD_MASK_SHAPEKEY | CD_MASK_RECAST |
-    CD_MASK_ORIGINDEX | CD_MASK_MVERT_SKIN | CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE;
+    CD_MASK_ORIGINDEX | CD_MASK_MVERT_SKIN | CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE |
+    CD_MASK_CUSTOMLOOPNORMAL;
 const CustomDataMask CD_MASK_BMESH =
     CD_MASK_MLOOPUV | CD_MASK_MLOOPCOL | CD_MASK_MTEXPOLY |
     CD_MASK_MSTICKY | CD_MASK_MDEFORMVERT | CD_MASK_PROP_FLT | CD_MASK_PROP_INT |
     CD_MASK_PROP_STR | CD_MASK_SHAPEKEY | CD_MASK_SHAPE_KEYINDEX | CD_MASK_MDISPS |
     CD_MASK_CREASE | CD_MASK_BWEIGHT | CD_MASK_RECAST | CD_MASK_PAINT_MASK |
-    CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN | CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE;
+    CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN | CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE |
+    CD_MASK_CUSTOMLOOPNORMAL;
 const CustomDataMask CD_MASK_FACECORNERS =  /* XXX Not used anywhere! */
     CD_MASK_MTFACE | CD_MASK_MCOL | CD_MASK_MTEXPOLY | CD_MASK_MLOOPUV |
     CD_MASK_MLOOPCOL | CD_MASK_NORMAL | CD_MASK_MLOOPTANGENT;
@@ -1316,7 +1322,7 @@ const CustomDataMask CD_MASK_EVERYTHING =
     /* BMESH ONLY END */
     CD_MASK_PAINT_MASK | CD_MASK_GRID_PAINT_MASK | CD_MASK_MVERT_SKIN |
     CD_MASK_FREESTYLE_EDGE | CD_MASK_FREESTYLE_FACE |
-    CD_MASK_MLOOPTANGENT | CD_MASK_TESSLOOPNORMAL;
+    CD_MASK_MLOOPTANGENT | CD_MASK_TESSLOOPNORMAL | CD_MASK_CUSTOMLOOPNORMAL;
 
 static const LayerTypeInfo *layerType_getInfo(int type)
 {
diff --git a/source/blender/blenkernel/intern/editderivedmesh.c b/source/blender/blenkernel/intern/editderivedmesh.c
index 62d9009a131..082edb01efd 100644
--- a/source/blender/blenkernel/intern/editderivedmesh.c
+++ b/source/blender/blenkernel/intern/editderivedmesh.c
@@ -169,12 +169,25 @@ static void emDM_calcNormals(DerivedMesh *dm)
 	dm->dirty &= ~DM_DIRTY_NORMALS;
 }
 
+static void emDM_calcLoopNormalsSpaceArray(
+        DerivedMesh *dm, const bool use_split_normals, const float split_angle, MLoopNorSpaceArray *r_lnors_spacearr);
+
 static void emDM_calcLoopNormals(DerivedMesh *dm, const bool use_split_normals, const float split_angle)
 {
+	emDM_calcLoopNormalsSpaceArray(dm, use_split_normals, split_angle, NULL);
+}
+
+/* #define DEBUG_CLNORS */
+
+static void emDM_calcLoopNormalsSpaceArray(
+        DerivedMesh *dm, const bool use_split_normals, const float split_angle, MLoopNorSpaceArray *r_lnors_spacearr)
+{
 	EditDerivedBMesh *bmdm = (EditDerivedBMesh *)dm;
 	BMesh *bm = bmdm->em->bm;
 	const float (*vertexCos)[3], (*vertexNos)[3], (*polyNos)[3];
 	float (*loopNos)[3];
+	short (*clnors_data)[2];
+	int cd_loop_clnors_offset;
 
 	/* calculate loop normals from poly and vertex normals */
 	emDM_ensureVertNormals(bmdm);
@@ -191,7 +204,37 @@ static void emDM_calcLoopNormals(DerivedMesh *dm, const bool use_split_normals,
 		loopNos = dm->getLoopDataArray(dm, CD_NORMAL);
 	}
 
-	BM_loops_calc_normal_vcos(bm, vertexCos, vertexNos, polyNos, use_split_normals, split_angle, loopNos);
+	/* We can have both, give priority to dm's data, and fallback to bm's ones. */
+	clnors_data = dm->getLoopDataArray(dm, CD_CUSTOMLOOPNORMAL);
+	cd_loop_clnors_offset = clnors_data ? -1 : CustomData_get_offset(&bm->ldata, CD_CUSTOMLOOPNORMAL);
+
+	BM_loops_calc_normal_vcos(bm, vertexCos, vertexNos, polyNos, use_split_normals, split_angle, loopNos,
+	                          r_lnors_spacearr, clnors_data, cd_loop_clnors_offset);
+#ifdef DEBUG_CLNORS
+	if (r_lnors_spacearr) {
+		int i;
+		for (i = 0; i < numLoops; i++) {
+			if (r_lnors_spacearr->lspacearr[i]->ref_alpha != 0.0f) {
+				LinkNode *loops = r_lnors_spacearr->lspacearr[i]->loops;
+				printf("Loop %d uses lnor space %p:\n", i, r_lnors_spacearr->lspacearr[i]);
+				print_v3("\tfinal lnor:", loopNos[i]);
+				print_v3("\tauto lnor:", r_lnors_spacearr->lspacearr[i]->vec_lnor);
+				print_v3("\tref_vec:", r_lnors_spacearr->lspacearr[i]->vec_ref);
+				printf("\talpha: %f\n\tbeta: %f\n\tloops: %p\n", r_lnors_spacearr->lspacearr[i]->ref_alpha,
+				       r_lnors_spacearr->lspacearr[i]->ref_beta, r_lnors_spacearr->lspacearr[i]->loops);
+				printf("\t\t(shared with loops");
+				while (loops) {
+					printf(" %d", GET_INT_FROM_POINTER(loops->link));
+					loops = loops->next;
+				}
+				printf(")\n");
+			}
+			else {
+				printf("Loop %d has no lnor space\n", i);
+			}
+		}
+	}
+#endif
 }
 
 static void emDM_recalcTessellation(DerivedMesh *UNUSED(dm))
@@ -1764,6 +1807,7 @@ DerivedMesh *getEditDerivedBMesh(BMEditMesh *em,
 
 	bmdm->dm.calcNormals = emDM_calcNormals;
 	bmdm->dm.calcLoopNormals = emDM_calcLoopNormals;
+	bmdm->dm.calcLoopNormalsSpaceArray = emDM_calcLoopNormalsSpaceArray;
 	bmdm->dm.recalcTessellation = emDM_recalcTessellation;
 
 	bmdm->dm.foreachMappedVert = emDM_foreachMappedVert;
diff --git a/source/blender/blenkernel/intern/mesh.c b/source/blender/blenkernel/intern/mesh.c
index 8ad20df18db..7757babaca2 100644
--- a/source/blender/blenkernel/intern/mesh.c
+++ b/source/blender/blenkernel/intern/mesh.c
@@ -414,6 +414,16 @@ void BKE_mesh_update_customdata_pointers(Mesh *me, const bool do_ensure_tess_cd)
 	me->mloopuv = CustomData_get_layer(&me->ldata, CD_MLOOPUV);
 }
 
+bool BKE_mesh_has_custom_loop_normals(Mesh *me)
+{
+	if (me->edit_btmesh) {
+		return CustomData_has_layer(&me->edit_btmesh->bm->ldata, CD_CUSTOMLOOPNORMAL);
+	}
+	else {
+		return CustomData_has_layer(&me->ldata, CD_CUSTOMLOOPNORMAL);
+	}
+}
+
 /* Note: unlinking is called when me->id.us is 0, question remains how
  * much unlinking of Library data in Mesh should be done... probably
  * we need a more generic method, like the expand() functions in
diff --git a/source/blender/blenkernel/intern/mesh_evaluate.c b/source/blender/blenkernel/intern/mesh_evaluate.c
index af9ac911921..f12fdeb7a80 100644
--- a/source/blender/blenkernel/intern/mesh_evaluate.c
+++ b/source/blender/blenkernel/intern/mesh_evaluate.c
@@ -39,6 +39,7 @@
 
 #include "BLI_utildefines.h"
 #include "BLI_memarena.h"
+#include "BLI_mempool.h"
 #include "BLI_math.h"
 #include "BLI_edgehash.h"
 #include "BLI_bitmap.h"
@@ -46,6 +47,8 @@
 #include "BLI_linklist.h"
 #include "BLI_linklist_stack.h"
 #include "BLI_alloca.h"
+#include "BLI_stack.h"
+#include "BLI_task.h"
 
 #include "BKE_customdata.h"
 #include "BKE_mesh.h"
@@ -58,8 +61,8 @@
 
 // #define DEBUG_TIME
 
+#include "PIL_time.h"
 #ifdef DEBUG_TIME
-#  include "PIL_time.h"
 #  include "PIL_time_utildefines.h"
 #endif
 
@@ -316,16 +319,751 @@ void BKE_mesh_calc_normals_tessface(MVert *mverts, int numVerts, MFace *mfaces,
 		MEM_freeN(fnors);
 }
 
+void BKE_lnor_spacearr_init(MLoopNorSpaceArray *lnors_spacearr, const int numLoops)
+{
+	if (!(lnors_spacearr->lspacearr && lnors_spacearr->loops_pool)) {
+		MemArena *mem;
+
+		if (!lnors_spacearr->mem) {
+			lnors_spacearr->mem = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
+		}
+		mem = lnors_spacearr->mem;
+		lnors_spacearr->lspacearr = BLI_memarena_calloc(mem, sizeof(MLoopNorSpace *) * (size_t)numLoops);
+		lnors_spacearr->loops_pool = BLI_memarena_alloc(mem, sizeof(LinkNode) * (size_t)numLoops);
+	}
+}
+
+void BKE_lnor_spacearr_clear(MLoopNorSpaceArray *lnors_spacearr)
+{
+	BLI_memarena_clear(lnors_spacearr->mem);
+	lnors_spacearr->lspacearr = NULL;
+	lnors_spacearr->loops_pool = NULL;
+}
+
+void BKE_lnor_spacearr_free(MLoopNorSpaceArray *lnors_spacearr)
+{
+	BLI_memarena_free(lnors_spacearr->mem);
+	lnors_spacearr->lspacearr = NULL;
+	lnors_spacearr->loops_pool = NULL;
+	lnors_spacearr->mem = NULL;
+}
+
+MLoopNorSpace *BKE_lnor_space_create(MLoopNorSpaceArray *lnors_spacearr)
+{
+	return BLI_memarena_calloc(lnors_spacearr->mem, sizeof(MLoopNorSpace));
+}
+
+/* This threshold is a bit touchy (usual float precision issue), this value seems OK. */
+#define LNOR_SPACE_TRIGO_THRESHOLD (1.0f - 1e-6f)
+
+/* Should only be called once.
+ * Beware, this modifies ref_vec and other_vec in place!
+ * In case no valid space can be generated, ref_alpha and ref_beta are set to zero (which means 'use auto lnors').
+ */
+void BKE_lnor_space_define(MLoopNorSpace *lnor_space, const float lnor[3],
+                           float vec_ref[3], float vec_other[3], BLI_Stack *edge_vectors)
+{
+	const float pi2 = (float)M_PI * 2.0f;
+	float tvec[3], dtp;
+	const float dtp_ref = dot_v3v3(vec_ref, lnor);
+	const float dtp_other = dot_v3v3(vec_other, lnor);
+
+	if (UNLIKELY(fabsf(dtp_ref) >= LNOR_SPACE_TRIGO_THRESHOLD || fabsf(dtp_other) >= LNOR_SPACE_TRIGO_THRESHOLD)) {
+		/* If vec_ref or vec_other are too much aligned with lnor, we can't build lnor space,
+		 * tag it as invalid and abort. */
+		lnor_space->ref_alpha = lnor_space->ref_beta = 0.0f;
+		return;
+	}
+
+	copy_v3_v3(lnor_space->vec_lnor, lnor);
+
+	/* Compute ref alpha, average angle of all available edge vectors to lnor. */
+	if (edge_vectors) {
+		float alpha = 0.0f;
+		int nbr = 0;
+		while (!BLI_stack_is_empty(edge_vectors)) {
+			const float *vec = BLI_stack_peek(edge_vectors);
+			alpha += saacosf(dot_v3v3(vec, lnor));
+			BLI_stack_discard(edge_vectors);
+			nbr++;
+		}
+		lnor_space->ref_alpha = alpha / (float)nbr;
+	}
+	else {
+		lnor_space->ref_alpha = (saacosf(dot_v3v3(vec_ref, lnor)) + saacosf(dot_v3v3(vec_other, lnor))) / 2.0f;
+	}
+
+	/* Project vec_ref on lnor's ortho plane. */
+	mul_v3_v3fl(tvec, lnor, dtp_ref);
+	sub_v3_v3(vec_ref, tvec);
+	normalize_v3_v3(lnor_space->vec_ref, vec_ref);
+
+	cross_v3_v3v3(tvec, lnor, lnor_space->vec_ref);
+	normalize_v3_v3(lnor_space->vec_ortho, tvec);
+
+	/* Project vec_other on lnor's ortho plane. */
+	mul_v3_v3fl(tvec, lnor, dtp_other);
+	sub_v3_v3(vec_other, tvec);
+	normalize_v3(vec_other);
+
+	/* Beta is angle between ref_vec and other_vec, around lnor. */
+	dtp = dot_v3v3(lnor_space->vec_ref, vec_other);
+	if (LIKELY(dtp < LNOR_SPACE_TRIGO_THRESHOLD)) {
+		const float beta = saacos(dtp);
+		lnor_space->ref_beta = (dot_v3v3(lnor_space->vec_ortho, vec_other) < 0.0f) ? pi2 - beta : beta;
+	}
+	else {
+		lnor_space->ref_beta = pi2;
+	}
+}
+
+void BKE_lnor_space_add_loop(MLoopNorSpaceArray *lnors_spacearr, MLoopNorSpace *lnor_space, const int ml_index,
+                             const bool do_add_loop)
+{
+	lnors_spacearr->lspacearr[ml_index] = lnor_space;
+	if (do_add_loop) {
+		BLI_linklist_prepend_nlink(&lnor_space->loops, SET_INT_IN_POINTER(ml_index), &lnors_spacearr->loops_pool[ml_index]);
+	}
+}
+
+MINLINE float unit_short_to_float(const short val)
+{
+	return (float)val / (float)SHRT_MAX;
+}
+
+MINLINE short unit_float_to_short(const float val)
+{
+	/* Rounding... */
+	return (short)floorf(val * (float)SHRT_MAX + 0.5f);
+}
+
+void BKE_lnor_space_custom_data_to_normal(MLoopNorSpace *lnor_space, const short clnor_data[2], float r_custom_lnor[3])
+{
+	/* NOP custom normal data or invalid lnor space, return. */
+	if (clnor_data[0] == 0 || lnor_space->ref_alpha == 0.0f || lnor_space->ref_beta == 0.0f) {
+		copy_v3_v3(r_custom_lnor, lnor_space->vec_lnor);
+		return;
+	}
+
+	{
+		/* TODO Check whether using sincosf() gives any noticeable benefit
+		 *      (could not even get it working under linux though)! */
+		const float pi2 = (float)(M_PI * 2.0);
+		const float alphafac = unit_short_to_float(clnor_data[0]);
+		const float alpha = (alphafac > 0.0f ? lnor_space->ref_alpha : pi2 - lnor_space->ref_alpha) * alphafac;
+		const float betafac = unit_short_to_float(clnor_data[1]);
+
+		mul_v3_v3fl(r_custom_lnor, lnor_space->vec_lnor, cosf(alpha));
+
+		if (betafac == 0.0f) {
+			madd_v3_v3fl(r_custom_lnor, lnor_space->vec_ref, sinf(alpha));
+		}
+		else {
+			const float sinalpha = sinf(alpha);
+			const float beta = (betafac > 0.0f ? lnor_space->ref_beta : pi2 - lnor_space->ref_beta) * betafac;
+			madd_v3_v3fl(r_custom_lnor, lnor_space->vec_ref, sinalpha * cosf(beta));
+			madd_v3_v3fl(r_custom_lnor, lnor_space->vec_ortho, sinalpha * sinf(beta));
+		}
+	}
+}
+
+void BKE_lnor_space_custom_normal_to_data(MLoopNorSpace *lnor_space, const float custom_lnor[3], short r_clnor_data[2])
+{
+	/* We use null vector as NOP custom normal (can be simpler than giving autocomputed lnor...). */
+	if (is_zero_v3(custom_lnor) || compare_v3v3(lnor_space->vec_lnor, custom_lnor, 1e-4f)) {
+		r_clnor_data[0] = r_clnor_data[1] = 0;
+		return;
+	}
+
+	{
+		const float pi2 = (float)(M_PI * 2.0);
+		const float cos_alpha = dot_v3v3(lnor_space->vec_lnor, custom_lnor);
+		float vec[3], cos_beta;
+		float alpha;
+
+		alpha = saacosf(cos_alpha);
+		if (alpha > lnor_space->ref_alpha) {
+			/* Note we could stick to [0, pi] range here, but makes decoding more complex, not worth it. */
+			r_clnor_data[0] = unit_float_to_short(-(pi2 - alpha) / (pi2 - lnor_space->ref_alpha));
+		}
+		else {
+			r_clnor_data[0] = unit_float_to_short(alpha / lnor_space->ref_alpha);
+		}
+
+		/* Project custom lnor on (vec_ref, vec_ortho) plane. */
+		mul_v3_v3fl(vec, lnor_space->vec_lnor, -cos_alpha);
+		add_v3_v3(vec, custom_lnor);
+		normalize_v3(vec);
+
+		cos_beta = dot_v3v3(lnor_space->vec_ref, vec);
+
+		if (cos_beta < LNOR_SPACE_TRIGO_THRESHOLD) {
+			float beta = saacosf(cos_beta);
+			if (dot_v3v3(lnor_space->vec_ortho, vec) < 0.0f) {
+				beta = pi2 - beta;
+			}
+
+			if (beta > lnor_space->ref_beta) {
+				r_clnor_data[1] = unit_float_to_short(-(pi2 - beta) / (pi2 - lnor_space->ref_beta));
+			}
+			else {
+				r_clnor_data[1] = unit_float_to_short(beta / lnor_space->ref_beta);
+			}
+		}
+		else {
+			r_clnor_data[1] = 0;
+		}
+	}
+}
+
+#define LOOP_SPLIT_TASK_BLOCK_SIZE 1024
+
+typedef struct LoopSplitTaskData {
+	/* Specific to each instance (each task). */
+	MLoopNorSpace *lnor_space;  /* We have to create those outside of tasks, since afaik memarena is not threadsafe. */
+	float (*lnor)[3];
+	const MLoop *ml_curr;
+	const MLoop *ml_prev;
+	int ml_curr_index;
+	int ml_prev_index;
+	const int *e2l_prev;  /* Also used a flag to switch between single or fan process! */
+	int mp_index;
+
+	/* This one is special, it's owned and managed by worker tasks, avoid to have to create it for each fan! */
+	BLI_Stack *edge_vectors;
+
+	char pad_c;
+} LoopSplitTaskData;
+
+typedef struct LoopSplitTaskDataCommon {
+	/* Read/write.
+	 * Note we do not need to protect it, though, since two different tasks will *always* affect different
+	 * elements in the arrays. */
+	MLoopNorSpaceArray *lnors_spacearr;
+	BLI_bitmap *sharp_verts;
+	float (*loopnors)[3];
+	short (*clnors_data)[2];
+
+	/* Read-only. */
+	const MVert *mverts;
+	const MEdge *medges;
+	const MLoop *mloops;
+	const MPoly *mpolys;
+	const int (*edge_to_loops)[2];
+	const int *loop_to_poly;
+	const float (*polynors)[3];
+
+	int numPolys;
+
+	/* ***** Workers communication. ***** */
+	ThreadQueue *task_queue;
+
+} LoopSplitTaskDataCommon;
+
+#define INDEX_UNSET INT_MIN
+#define INDEX_INVALID -1
+/* See comment about edge_to_loops below. */
+#define IS_EDGE_SHARP(_e2l) (ELEM((_e2l)[1], INDEX_UNSET, INDEX_INVALID))
+
+static void split_loop_nor_single_do(LoopSplitTaskDataCommon *common_data, LoopSplitTaskData *data)
+{
+	MLoopNorSpaceArray *lnors_spacearr = common_data->lnors_spacearr;
+	short (*clnors_data)[2] = common_data->clnors_data;
+
+	const MVert *mverts = common_data->mverts;
+	const MEdge *medges = common_data->medges;
+	const float (*polynors)[3] = common_data->polynors;
+
+	MLoopNorSpace *lnor_space = data->lnor_space;
+	float (*lnor)[3] = data->lnor;
+	const MLoop *ml_curr = data->ml_curr;
+	const MLoop *ml_prev = data->ml_prev;
+	const int ml_curr_index = data->ml_curr_index;
+#if 0  /* Not needed for 'single' loop. */
+	const int ml_prev_index = data->ml_prev_index;
+	const int *e2l_prev = data->e2l_prev;
+#endif
+	const int mp_index = data->mp_index;
+
+	/* Simple case (both edges around that vertex are sharp in current polygon),
+	 * this loop just takes its poly normal.
+	 */
+	copy_v3_v3(*lnor, polynors[mp_index]);
+
+	/* printf("BASIC: handling loop %d / edge %d / vert %d\n", ml_curr_index, ml_curr->e, ml_curr->v); */
+
+	/* If needed, generate this (simple!) lnor space. */
+	if (lnors_spacearr) {
+		float vec_curr[3], vec_prev[3];
+
+		const unsigned int mv_pivot_index = ml_curr->v;  /* The vertex we are "fanning" around! */
+		const MVert *mv_pivot = &mverts[mv_pivot_index];
+		const MEdge *me_curr = &medges[ml_curr->e];
+		const MVert *mv_2 = (me_curr->v1 == mv_pivot_index) ? &mverts[me_curr->v2] : &mverts[me_curr->v1];
+		const MEdge *me_prev = &medges[ml_prev->e];
+		const MVert *mv_3 = (me_prev->v1 == mv_pivot_index) ? &mverts[me_prev->v2] : &mverts[me_prev->v1];
+
+		sub_v3_v3v3(vec_curr, mv_2->co, mv_pivot->co);
+		normalize_v3(vec_curr);
+		sub_v3_v3v3(vec_prev, mv_3->co, mv_pivot->co);
+		normalize_v3(vec_prev);
+
+		BKE_lnor_space_define(lnor_space, *lnor, vec_curr, vec_prev, NULL);
+		/* We know there is only one loop in this space, no need to create a linklist in this case... */
+		BKE_lnor_space_add_loop(lnors_spacearr, lnor_space, ml_curr_index, false);
+
+		if (clnors_data) {
+			BKE_lnor_space_custom_data_to_normal(lnor_space, clnors_data[ml_curr_index], *lnor);
+		}
+	}
+}
+
+static void split_loop_nor_fan_do(LoopSplitTaskDataCommon *common_data, LoopSplitTaskData *data)
+{
+	MLoopNorSpaceArray *lnors_spacearr = common_data->lnors_spacearr;
+	float (*loopnors)[3] = common_data->loopnors;
+	short (*clnors_data)[2] = common_data->clnors_data;
+
+	const MVert *mverts = common_data->mverts;
+	const MEdge *medges = common_data->medges;
+	const MLoop *mloops = common_data->mloops;
+	const MPoly *mpolys = common_data->mpolys;
+	const int (*edge_to_loops)[2] = common_data->edge_to_loops;
+	const int *loop_to_poly = common_data->loop_to_poly;
+	const float (*polynors)[3] = common_data->polynors;
+
+	MLoopNorSpace *lnor_space = data->lnor_space;
+#if 0  /* Not needed for 'fan' loops. */
+	float (*lnor)[3] = data->lnor;
+#endif
+	const MLoop *ml_curr = data->ml_curr;
+	const MLoop *ml_prev = data->ml_prev;
+	const int ml_curr_index = data->ml_curr_index;
+	const int ml_prev_index = data->ml_prev_index;
+	const int mp_index = data->mp_index;
+	const int *e2l_prev = data->e2l_prev;
+
+	BLI_Stack *edge_vectors = data->edge_vectors;
+
+	/* Gah... We have to fan around current vertex, until we find the other non-smooth edge,
+	 * and accumulate face normals into the vertex!
+	 * Note in case this vertex has only one sharp edges, this is a waste because the normal is the same as
+	 * the vertex normal, but I do not see any easy way to detect that (would need to count number
+	 * of sharp edges per vertex, I doubt the additional memory usage would be worth it, especially as
+	 * it should not be a common case in real-life meshes anyway).
+	 */
+	const unsigned int mv_pivot_index = ml_curr->v;  /* The vertex we are "fanning" around! */
+	const MVert *mv_pivot = &mverts[mv_pivot_index];
+	const MEdge *me_org = &medges[ml_curr->e];  /* ml_curr would be mlfan_prev if we needed that one */
+	const int *e2lfan_curr;
+	float vec_curr[3], vec_prev[3], vec_org[3];
+	const MLoop *mlfan_curr, *mlfan_next;
+	const MPoly *mpfan_next;
+	float lnor[3] = {0.0f, 0.0f, 0.0f};
+	/* mlfan_vert_index: the loop of our current edge might not be the loop of our current vertex! */
+	int mlfan_curr_index, mlfan_vert_index, mpfan_curr_index;
+
+	/* We validate clnors data on the fly - cheapest way to do! */
+	int clnors_avg[2] = {0, 0};
+	short (*clnor_ref)[2] = NULL;
+	int clnors_nbr = 0;
+	bool clnors_invalid = false;
+
+	/* Temp loop normal stack. */
+	BLI_SMALLSTACK_DECLARE(normal, float *);
+	/* Temp clnors stack. */
+	BLI_SMALLSTACK_DECLARE(clnors, short *);
+
+	e2lfan_curr = e2l_prev;
+	mlfan_curr = ml_prev;
+	mlfan_curr_index = ml_prev_index;
+	mlfan_vert_index = ml_curr_index;
+	mpfan_curr_index = mp_index;
+
+	BLI_assert(mlfan_curr_index >= 0);
+	BLI_assert(mlfan_vert_index >= 0);
+	BLI_assert(mpfan_curr_index >= 0);
+
+	/* Only need to compute previous edge's vector once, then we can just reuse old current one! */
+	{
+		const MVert *mv_2 = (me_org->v1 == mv_pivot_index) ? &mverts[me_org->v2] : &mverts[me_org->v1];
+
+		sub_v3_v3v3(vec_org, mv_2->co, mv_pivot->co);
+		normalize_v3(vec_org);
+		copy_v3_v3(vec_prev, vec_org);
+
+		if (lnors_spacearr) {
+			BLI_stack_push(edge_vectors, vec_org);
+		}
+	}
+
+	/* printf("FAN: vert %d, start edge %d\n", mv_pivot_index, ml_curr->e); */
+
+	while (true) {
+		const MEdge *me_curr = &medges[mlfan_curr->e];
+		/* Compute edge vectors.
+		 * NOTE: We could pre-compute those into an array, in the first iteration, instead of computing them
+		 *       twice (or more) here. However, time gained is not worth memory and time lost,
+		 *       given the fact that this code should not be called that much in real-life meshes...
+		 */
+		{
+			const MVert *mv_2 = (me_curr->v1 == mv_pivot_index) ? &mverts[me_curr->v2] : &mverts[me_curr->v1];
+
+			sub_v3_v3v3(vec_curr, mv_2->co, mv_pivot->co);
+			normalize_v3(vec_curr);
+		}
+
+		/* printf("\thandling edge %d / loop %d\n", mlfan_curr->e, mlfan_curr_index); */
+
+		{
+			/* Code similar to accumulate_vertex_normals_poly. */
+			/* Calculate angle between the two poly edges incident on this vertex. */
+			const float fac = saacos(dot_v3v3(vec_curr, vec_prev));
+			/* Accumulate */
+			madd_v3_v3fl(lnor, polynors[mpfan_curr_index], fac);
+
+			if (clnors_data) {
+				/* Accumulate all clnors, if they are not all equal we have to fix that! */
+				short (*clnor)[2] = &clnors_data[mlfan_vert_index];
+				if (clnors_nbr) {
+					clnors_invalid |= ((*clnor_ref)[0] != (*clnor)[0] || (*clnor_ref)[1] != (*clnor)[1]);
+				}
+				else {
+					clnor_ref = clnor;
+				}
+				clnors_avg[0] += (*clnor)[0];
+				clnors_avg[1] += (*clnor)[1];
+				clnors_nbr++;
+				/* We store here a pointer to all custom lnors processed. */
+				BLI_SMALLSTACK_PUSH(clnors, (short *)*clnor);
+			}
+		}
+
+		/* We store here a pointer to all loop-normals processed. */
+		BLI_SMALLSTACK_PUSH(normal, (float *)(loopnors[mlfan_vert_index]));
+
+		if (lnors_spacearr) {
+			/* Assign current lnor space to current 'vertex' loop. */
+			BKE_lnor_space_add_loop(lnors_spacearr, lnor_space, mlfan_vert_index, true);
+			if (me_curr != me_org) {
+				/* We store here all edges-normalized vectors processed. */
+				BLI_stack_push(edge_vectors, vec_curr);
+			}
+		}
+
+		if (IS_EDGE_SHARP(e2lfan_curr) || (me_curr == me_org)) {
+			/* Current edge is sharp and we have finished with this fan of faces around this vert,
+			 * or this vert is smooth, and we have completed a full turn around it.
+			 */
+			/* printf("FAN: Finished!\n"); */
+			break;
+		}
+
+		copy_v3_v3(vec_prev, vec_curr);
+
+		/* Warning! This is rather complex!
+		 * We have to find our next edge around the vertex (fan mode).
+		 * First we find the next loop, which is either previous or next to mlfan_curr_index, depending
+		 * whether both loops using current edge are in the same direction or not, and whether
+		 * mlfan_curr_index actually uses the vertex we are fanning around!
+		 * mlfan_curr_index is the index of mlfan_next here, and mlfan_next is not the real next one
+		 * (i.e. not the future mlfan_curr)...
+		 */
+		mlfan_curr_index = (e2lfan_curr[0] == mlfan_curr_index) ? e2lfan_curr[1] : e2lfan_curr[0];
+		mpfan_curr_index = loop_to_poly[mlfan_curr_index];
+
+		BLI_assert(mlfan_curr_index >= 0);
+		BLI_assert(mpfan_curr_index >= 0);
+
+		mlfan_next = &mloops[mlfan_curr_index];
+		mpfan_next = &mpolys[mpfan_curr_index];
+		if ((mlfan_curr->v == mlfan_next->v && mlfan_curr->v == mv_pivot_index) ||
+		    (mlfan_curr->v != mlfan_next->v && mlfan_curr->v != mv_pivot_index))
+		{
+			/* We need the previous loop, but current one is our vertex's loop. */
+			mlfan_vert_index = mlfan_curr_index;
+			if (--mlfan_curr_index < mpfan_next->loopstart) {
+				mlfan_curr_index = mpfan_next->loopstart + mpfan_next->totloop - 1;
+			}
+		}
+		else {
+			/* We need the next loop, which is also our vertex's loop. */
+			if (++mlfan_curr_index >= mpfan_next->loopstart + mpfan_next->totloop) {
+				mlfan_curr_index = mpfan_next->loopstart;
+			}
+			mlfan_vert_index = mlfan_curr_index;
+		}
+		mlfan_curr = &mloops[mlfan_curr_index];
+		/* And now we are back in sync, mlfan_curr_index is the index of mlfan_curr! Pff! */
+
+		e2lfan_curr = edge_to_loops[mlfan_curr->e];
+	}
+
+	{
+		float lnor_len = normalize_v3(lnor);
+
+		/* If we are generating lnor spacearr, we can now define the one for this fan,
+		 * and optionally compute final lnor from custom data too!
+		 */
+		if (lnors_spacearr) {
+			if (UNLIKELY(lnor_len == 0.0f)) {
+				/* Use vertex normal as fallback! */
+				copy_v3_v3(lnor, loopnors[mlfan_vert_index]);
+				lnor_len = 1.0f;
+			}
+
+			BKE_lnor_space_define(lnor_space, lnor, vec_org, vec_curr, edge_vectors);
+
+			if (clnors_data) {
+				if (clnors_invalid) {
+					short *clnor;
+
+					clnors_avg[0] /= clnors_nbr;
+					clnors_avg[1] /= clnors_nbr;
+					/* Fix/update all clnors of this fan with computed average value. */
+					printf("Invalid clnors in this fan!\n");
+					while ((clnor = BLI_SMALLSTACK_POP(clnors))) {
+						//print_v2("org clnor", clnor);
+						clnor[0] = (short)clnors_avg[0];
+						clnor[1] = (short)clnors_avg[1];
+					}
+					//print_v2("new clnors", clnors_avg);
+				}
+				/* Extra bonus: since smallstack is local to this func, no more need to empty it at all cost! */
+
+				BKE_lnor_space_custom_data_to_normal(lnor_space, *clnor_ref, lnor);
+			}
+		}
+
+		/* In case we get a zero normal here, just use vertex normal already set! */
+		if (LIKELY(lnor_len != 0.0f)) {
+			/* Copy back the final computed normal into all related loop-normals. */
+			float *nor;
+
+			while ((nor = BLI_SMALLSTACK_POP(normal))) {
+				copy_v3_v3(nor, lnor);
+			}
+		}
+		/* Extra bonus: since smallstack is local to this func, no more need to empty it at all cost! */
+	}
+}
+
+static void loop_split_worker_do(
+        LoopSplitTaskDataCommon *common_data, LoopSplitTaskData *data, BLI_Stack *edge_vectors)
+{
+	BLI_assert(data->ml_curr);
+	if (data->e2l_prev) {
+		BLI_assert((edge_vectors == NULL) || BLI_stack_is_empty(edge_vectors));
+		data->edge_vectors = edge_vectors;
+		split_loop_nor_fan_do(common_data, data);
+	}
+	else {
+		/* No need for edge_vectors for 'single' case! */
+		split_loop_nor_single_do(common_data, data);
+	}
+}
+
+static void loop_split_worker(TaskPool *UNUSED(pool), void *taskdata, int UNUSED(threadid))
+{
+	LoopSplitTaskDataCommon *common_data = taskdata;
+	LoopSplitTaskData *data_buff;
+
+	/* Temp edge vectors stack, only used when computing lnor spacearr. */
+	BLI_Stack *edge_vectors = common_data->lnors_spacearr ? BLI_stack_new(sizeof(float[3]), __func__) : NULL;
+
+#ifdef DEBUG_TIME
+	TIMEIT_START(loop_split_worker);
+#endif
+
+	while ((data_buff = BLI_thread_queue_pop(common_data->task_queue))) {
+		LoopSplitTaskData *data = data_buff;
+		int i;
+
+		for (i = 0; i < LOOP_SPLIT_TASK_BLOCK_SIZE; i++, data++) {
+			/* A NULL ml_curr is used to tag ended data! */
+			if (data->ml_curr == NULL) {
+				break;
+			}
+			loop_split_worker_do(common_data, data, edge_vectors);
+		}
+
+		MEM_freeN(data_buff);
+	}
+
+	if (edge_vectors) {
+		BLI_stack_free(edge_vectors);
+	}
+
+#ifdef DEBUG_TIME
+	TIMEIT_END(loop_split_worker);
+#endif
+}
+
+/* Note we use data_buff to detect whether we are in threaded context or not, in later case it is NULL. */
+static void loop_split_generator_do(LoopSplitTaskDataCommon *common_data, const bool threaded)
+{
+	MLoopNorSpaceArray *lnors_spacearr = common_data->lnors_spacearr;
+	BLI_bitmap *sharp_verts = common_data->sharp_verts;
+	float (*loopnors)[3] = common_data->loopnors;
+
+	const MLoop *mloops = common_data->mloops;
+	const MPoly *mpolys = common_data->mpolys;
+	const int (*edge_to_loops)[2] = common_data->edge_to_loops;
+	const int numPolys = common_data->numPolys;
+
+	const MPoly *mp;
+	int mp_index;
+
+	LoopSplitTaskData *data, *data_buff = NULL, data_mem;
+	int data_idx = 0;
+
+	/* Temp edge vectors stack, only used when computing lnor spacearr (and we are not multi-threading). */
+	BLI_Stack *edge_vectors = (lnors_spacearr && !data_buff) ? BLI_stack_new(sizeof(float[3]), __func__) : NULL;
+
+#ifdef DEBUG_TIME
+	TIMEIT_START(loop_split_generator);
+#endif
+
+	if (!threaded) {
+		memset(&data_mem, 0, sizeof(data_mem));
+		data = &data_mem;
+	}
+
+	/* We now know edges that can be smoothed (with their vector, and their two loops), and edges that will be hard!
+	 * Now, time to generate the normals.
+	 */
+	for (mp = mpolys, mp_index = 0; mp_index < numPolys; mp++, mp_index++) {
+		const MLoop *ml_curr, *ml_prev;
+		float (*lnors)[3];
+		const int ml_last_index = (mp->loopstart + mp->totloop) - 1;
+		int ml_curr_index = mp->loopstart;
+		int ml_prev_index = ml_last_index;
+
+		ml_curr = &mloops[ml_curr_index];
+		ml_prev = &mloops[ml_prev_index];
+		lnors = &loopnors[ml_curr_index];
+
+		for (; ml_curr_index <= ml_last_index; ml_curr++, ml_curr_index++, lnors++) {
+			const int *e2l_curr = edge_to_loops[ml_curr->e];
+			const int *e2l_prev = edge_to_loops[ml_prev->e];
+
+			if (!IS_EDGE_SHARP(e2l_curr) && (!lnors_spacearr || BLI_BITMAP_TEST_BOOL(sharp_verts, ml_curr->v))) {
+				/* A smooth edge, and we are not generating lnor_spacearr, or the related vertex is sharp.
+				 * We skip it because it is either:
+				 * - in the middle of a 'smooth fan' already computed (or that will be as soon as we hit
+				 *   one of its ends, i.e. one of its two sharp edges), or...
+				 * - the related vertex is a "full smooth" one, in which case pre-populated normals from vertex
+				 *   are just fine (or it has already be handled in a previous loop in case of needed lnors spacearr)!
+				 */
+				/* printf("Skipping loop %d / edge %d / vert %d(%d)\n", ml_curr_index, ml_curr->e, ml_curr->v, sharp_verts[ml_curr->v]); */
+			}
+			else {
+				if (threaded) {
+					if (data_idx == 0) {
+						data_buff = MEM_callocN(sizeof(*data_buff) * LOOP_SPLIT_TASK_BLOCK_SIZE, __func__);
+					}
+					data = &data_buff[data_idx];
+				}
+
+				if (IS_EDGE_SHARP(e2l_curr) && IS_EDGE_SHARP(e2l_prev)) {
+					data->lnor = lnors;
+					data->ml_curr = ml_curr;
+					data->ml_prev = ml_prev;
+					data->ml_curr_index = ml_curr_index;
+#if 0  /* Not needed for 'single' loop. */
+					data->ml_prev_index = ml_prev_index;
+					data->e2l_prev = NULL;  /* Tag as 'single' task. */
+#endif
+					data->mp_index = mp_index;
+					if (lnors_spacearr) {
+						data->lnor_space = BKE_lnor_space_create(lnors_spacearr);
+					}
+				}
+				/* We *do not need* to check/tag loops as already computed!
+				 * Due to the fact a loop only links to one of its two edges, a same fan *will never be walked
+				 * more than once!*
+				 * Since we consider edges having neighbor polys with inverted (flipped) normals as sharp, we are sure
+				 * that no fan will be skipped, even only considering the case (sharp curr_edge, smooth prev_edge),
+				 * and not the alternative (smooth curr_edge, sharp prev_edge).
+				 * All this due/thanks to link between normals and loop ordering (i.e. winding).
+				 */
+				else {
+#if 0  /* Not needed for 'fan' loops. */
+					data->lnor = lnors;
+#endif
+					data->ml_curr = ml_curr;
+					data->ml_prev = ml_prev;
+					data->ml_curr_index = ml_curr_index;
+					data->ml_prev_index = ml_prev_index;
+					data->e2l_prev = e2l_prev;  /* Also tag as 'fan' task. */
+					data->mp_index = mp_index;
+					if (lnors_spacearr) {
+						data->lnor_space = BKE_lnor_space_create(lnors_spacearr);
+						/* Tag related vertex as sharp, to avoid fanning around it again (in case it was a smooth one).
+						 * This *has* to be done outside of workers tasks! */
+						BLI_BITMAP_ENABLE(sharp_verts, ml_curr->v);
+					}
+				}
+
+				if (threaded) {
+					data_idx++;
+					if (data_idx == LOOP_SPLIT_TASK_BLOCK_SIZE) {
+						BLI_thread_queue_push(common_data->task_queue, data_buff);
+						data_idx = 0;
+					}
+				}
+				else {
+					loop_split_worker_do(common_data, data, edge_vectors);
+					memset(data, 0, sizeof(data_mem));
+				}
+			}
+
+			ml_prev = ml_curr;
+			ml_prev_index = ml_curr_index;
+		}
+	}
+
+	if (threaded) {
+		/* Last block of data... Since it is calloc'ed and we use first NULL item as stopper, everything is fine. */
+		if (LIKELY(data_idx)) {
+			BLI_thread_queue_push(common_data->task_queue, data_buff);
+		}
+
+		/* This will signal all other worker threads to wake up and finish! */
+		BLI_thread_queue_nowait(common_data->task_queue);
+	}
+
+	if (edge_vectors) {
+		BLI_stack_free(edge_vectors);
+	}
+
+#ifdef DEBUG_TIME
+	TIMEIT_END(loop_split_generator);
+#endif
+}
+
+static void loop_split_generator(TaskPool *UNUSED(pool), void *taskdata, int UNUSED(threadid))
+{
+	LoopSplitTaskDataCommon *common_data = taskdata;
+
+	loop_split_generator_do(common_data, true);
+}
+
 /**
  * Compute split normals, i.e. vertex normals associated with each poly (hence 'loop normals').
  * Useful to materialize sharp edges (or non-smooth faces) without actually modifying the geometry (splitting edges).
  */
 void BKE_mesh_normals_loop_split(
-        MVert *mverts, const int UNUSED(numVerts), MEdge *medges, const int numEdges,
+        MVert *mverts, const int numVerts, MEdge *medges, const int numEdges,
         MLoop *mloops, float (*r_loopnors)[3], const int numLoops,
-        MPoly *mpolys, float (*polynors)[3], const int numPolys,
-        const bool use_split_normals, float split_angle)
+        MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        const bool use_split_normals, float split_angle,
+        MLoopNorSpaceArray *r_lnors_spacearr, short (*clnors_data)[2], int *r_loop_to_poly)
 {
+
+	/* For now this is not supported. If we do not use split normals, we do not generate anything fancy! */
+	BLI_assert(use_split_normals || !(r_lnors_spacearr || r_loop_to_poly));
+
 	if (!use_split_normals) {
 		/* In this case, we simply fill lnors with vnors, quite simple!
 		 * Note this is done here to keep some logic and consistency in this quite complex code,
@@ -342,11 +1080,6 @@ void BKE_mesh_normals_loop_split(
 
 	{
 
-#define INDEX_UNSET INT_MIN
-#define INDEX_INVALID -1
-/* See comment about edge_to_loops below. */
-#define IS_EDGE_SHARP(_e2l) (ELEM((_e2l)[1], INDEX_UNSET, INDEX_INVALID))
-
 	/* Mapping edge -> loops.
 	 * If that edge is used by more than two loops (polys), it is always sharp (and tagged as such, see below).
 	 * We also use the second loop index as a kind of flag: smooth edge: > 0,
@@ -359,21 +1092,39 @@ void BKE_mesh_normals_loop_split(
 	int (*edge_to_loops)[2] = MEM_callocN(sizeof(int[2]) * (size_t)numEdges, __func__);
 
 	/* Simple mapping from a loop to its polygon index. */
-	int *loop_to_poly = MEM_mallocN(sizeof(int) * (size_t)numLoops, __func__);
+	int *loop_to_poly = r_loop_to_poly ? r_loop_to_poly : MEM_mallocN(sizeof(int) * (size_t)numLoops, __func__);
 
 	MPoly *mp;
-	int mp_index;
-	const bool check_angle = (split_angle < (float)M_PI);
+	int mp_index, me_index;
+	bool check_angle = (split_angle < (float)M_PI);
+	int i;
 
-	/* Temp normal stack. */
-	BLI_SMALLSTACK_DECLARE(normal, float *);
+	BLI_bitmap *sharp_verts = NULL;
+	MLoopNorSpaceArray _lnors_spacearr = {NULL};
+
+	LoopSplitTaskDataCommon common_data = {NULL};
 
 #ifdef DEBUG_TIME
 	TIMEIT_START(BKE_mesh_normals_loop_split);
 #endif
 
 	if (check_angle) {
-		split_angle = cosf(split_angle);
+		/* When using custom loop normals, disable the angle feature! */
+		if (clnors_data) {
+			check_angle = false;
+		}
+		else {
+			split_angle = cosf(split_angle);
+		}
+	}
+
+	if (!r_lnors_spacearr && clnors_data) {
+		/* We need to compute lnor spacearr if some custom lnor data are given to us! */
+		r_lnors_spacearr = &_lnors_spacearr;
+	}
+	if (r_lnors_spacearr) {
+		BKE_lnor_spacearr_init(r_lnors_spacearr, numLoops);
+		sharp_verts = BLI_BITMAP_NEW((size_t)numVerts, __func__);
 	}
 
 	/* This first loop check which edges are actually smooth, and compute edge vectors. */
@@ -427,189 +1178,271 @@ void BKE_mesh_normals_loop_split(
 		}
 	}
 
-	/* We now know edges that can be smoothed (with their vector, and their two loops), and edges that will be hard!
-	 * Now, time to generate the normals.
-	 */
-	for (mp = mpolys, mp_index = 0; mp_index < numPolys; mp++, mp_index++) {
-		MLoop *ml_curr, *ml_prev;
-		float (*lnors)[3];
-		const int ml_last_index = (mp->loopstart + mp->totloop) - 1;
-		int ml_curr_index = mp->loopstart;
-		int ml_prev_index = ml_last_index;
+	if (r_lnors_spacearr) {
+		/* Tag vertices that have at least one sharp edge as 'sharp' (used for the lnor spacearr computation).
+		 * XXX This third loop over edges is a bit disappointing, could not find any other way yet.
+		 *     Not really performance-critical anyway.
+		 */
+		for (me_index = 0; me_index < numEdges; me_index++) {
+			const int *e2l = edge_to_loops[me_index];
+			const MEdge *me = &medges[me_index];
+			if (IS_EDGE_SHARP(e2l)) {
+				BLI_BITMAP_ENABLE(sharp_verts, me->v1);
+				BLI_BITMAP_ENABLE(sharp_verts, me->v2);
+			}
+		}
+	}
 
-		ml_curr = &mloops[ml_curr_index];
-		ml_prev = &mloops[ml_prev_index];
-		lnors = &r_loopnors[ml_curr_index];
+	/* Init data common to all tasks. */
+	common_data.lnors_spacearr = r_lnors_spacearr;
+	common_data.loopnors = r_loopnors;
+	common_data.clnors_data = clnors_data;
 
-		for (; ml_curr_index <= ml_last_index; ml_curr++, ml_curr_index++, lnors++) {
-			const int *e2l_curr = edge_to_loops[ml_curr->e];
-			const int *e2l_prev = edge_to_loops[ml_prev->e];
+	common_data.mverts = mverts;
+	common_data.medges = medges;
+	common_data.mloops = mloops;
+	common_data.mpolys = mpolys;
+	common_data.sharp_verts = sharp_verts;
+	common_data.edge_to_loops = (const int(*)[2])edge_to_loops;
+	common_data.loop_to_poly = loop_to_poly;
+	common_data.polynors = polynors;
+	common_data.numPolys = numPolys;
 
-			if (!IS_EDGE_SHARP(e2l_curr)) {
-				/* A smooth edge.
-				 * We skip it because it is either:
-				 * - in the middle of a 'smooth fan' already computed (or that will be as soon as we hit
-				 *   one of its ends, i.e. one of its two sharp edges), or...
-				 * - the related vertex is a "full smooth" one, in which case pre-populated normals from vertex
-				 *   are just fine!
-				 */
-			}
-			else if (IS_EDGE_SHARP(e2l_prev)) {
-				/* Simple case (both edges around that vertex are sharp in current polygon),
-				 * this vertex just takes its poly normal.
-				 */
-				copy_v3_v3(*lnors, polynors[mp_index]);
-				/* No need to mark loop as done here, we won't run into it again anyway! */
-			}
-			/* We *do not need* to check/tag loops as already computed!
-			 * Due to the fact a loop only links to one of its two edges, a same fan *will never be walked more than
-			 * once!*
-			 * Since we consider edges having neighbor polys with inverted (flipped) normals as sharp, we are sure that
-			 * no fan will be skipped, even only considering the case (sharp curr_edge, smooth prev_edge), and not the
-			 * alternative (smooth curr_edge, sharp prev_edge).
-			 * All this due/thanks to link between normals and loop ordering.
-			 */
-			else {
-				/* Gah... We have to fan around current vertex, until we find the other non-smooth edge,
-				 * and accumulate face normals into the vertex!
-				 * Note in case this vertex has only one sharp edges, this is a waste because the normal is the same as
-				 * the vertex normal, but I do not see any easy way to detect that (would need to count number
-				 * of sharp edges per vertex, I doubt the additional memory usage would be worth it, especially as
-				 * it should not be a common case in real-life meshes anyway).
-				 */
-				const unsigned int mv_pivot_index = ml_curr->v;  /* The vertex we are "fanning" around! */
-				const MVert *mv_pivot = &mverts[mv_pivot_index];
-				const int *e2lfan_curr;
-				float vec_curr[3], vec_prev[3];
-				MLoop *mlfan_curr, *mlfan_next;
-				MPoly *mpfan_next;
-				float lnor[3] = {0.0f, 0.0f, 0.0f};
-				/* mlfan_vert_index: the loop of our current edge might not be the loop of our current vertex! */
-				int mlfan_curr_index, mlfan_vert_index, mpfan_curr_index;
-
-				e2lfan_curr = e2l_prev;
-				mlfan_curr = ml_prev;
-				mlfan_curr_index = ml_prev_index;
-				mlfan_vert_index = ml_curr_index;
-				mpfan_curr_index = mp_index;
-
-				BLI_assert(mlfan_curr_index >= 0);
-				BLI_assert(mlfan_vert_index >= 0);
-				BLI_assert(mpfan_curr_index >= 0);
-
-				/* Only need to compute previous edge's vector once, then we can just reuse old current one! */
-				{
-					const MEdge *me_prev = &medges[ml_curr->e];  /* ml_curr would be mlfan_prev if we needed that one */
-					const MVert *mv_2 = (me_prev->v1 == mv_pivot_index) ? &mverts[me_prev->v2] : &mverts[me_prev->v1];
+	if (numLoops < LOOP_SPLIT_TASK_BLOCK_SIZE * 8) {
+		/* Not enough loops to be worth the whole threading overhead... */
+		loop_split_generator_do(&common_data, false);
+	}
+	else {
+		TaskScheduler *task_scheduler;
+		TaskPool *task_pool;
+		int nbr_workers;
 
-					sub_v3_v3v3(vec_prev, mv_2->co, mv_pivot->co);
-					normalize_v3(vec_prev);
-				}
+		common_data.task_queue = BLI_thread_queue_init();
 
-				while (true) {
-					/* Compute edge vectors.
-					 * NOTE: We could pre-compute those into an array, in the first iteration, instead of computing them
-					 *       twice (or more) here. However, time gained is not worth memory and time lost,
-					 *       given the fact that this code should not be called that much in real-life meshes...
-					 */
-					{
-						const MEdge *me_curr = &medges[mlfan_curr->e];
-						const MVert *mv_2 = (me_curr->v1 == mv_pivot_index) ? &mverts[me_curr->v2] :
-						                                                      &mverts[me_curr->v1];
-
-						sub_v3_v3v3(vec_curr, mv_2->co, mv_pivot->co);
-						normalize_v3(vec_curr);
-					}
+		task_scheduler = BLI_task_scheduler_get();
+		task_pool = BLI_task_pool_create(task_scheduler, NULL);
 
-					{
-						/* Code similar to accumulate_vertex_normals_poly. */
-						/* Calculate angle between the two poly edges incident on this vertex. */
-						const float fac = saacos(dot_v3v3(vec_curr, vec_prev));
-						/* Accumulate */
-						madd_v3_v3fl(lnor, polynors[mpfan_curr_index], fac);
-					}
+		nbr_workers = max_ii(2, BLI_task_scheduler_num_threads(task_scheduler));
+		for (i = 1; i < nbr_workers; i++) {
+			BLI_task_pool_push(task_pool, loop_split_worker, &common_data, false, TASK_PRIORITY_HIGH);
+		}
+		BLI_task_pool_push(task_pool, loop_split_generator, &common_data, false, TASK_PRIORITY_HIGH);
+		BLI_task_pool_work_and_wait(task_pool);
 
-					/* We store here a pointer to all loop-normals processed. */
-					BLI_SMALLSTACK_PUSH(normal, &(r_loopnors[mlfan_vert_index][0]));
+		BLI_task_pool_free(task_pool);
 
-					if (IS_EDGE_SHARP(e2lfan_curr)) {
-						/* Current edge is sharp, we have finished with this fan of faces around this vert! */
-						break;
-					}
+		BLI_thread_queue_free(common_data.task_queue);
+	}
+
+	MEM_freeN(edge_to_loops);
+	if (!r_loop_to_poly) {
+		MEM_freeN(loop_to_poly);
+	}
+
+	if (r_lnors_spacearr) {
+		MEM_freeN(sharp_verts);
+		if (r_lnors_spacearr == &_lnors_spacearr) {
+			BKE_lnor_spacearr_free(r_lnors_spacearr);
+		}
+	}
+
+#ifdef DEBUG_TIME
+	TIMEIT_END(BKE_mesh_normals_loop_split);
+#endif
+
+	}
+}
 
-					copy_v3_v3(vec_prev, vec_curr);
-
-					/* Warning! This is rather complex!
-					 * We have to find our next edge around the vertex (fan mode).
-					 * First we find the next loop, which is either previous or next to mlfan_curr_index, depending
-					 * whether both loops using current edge are in the same direction or not, and whether
-					 * mlfan_curr_index actually uses the vertex we are fanning around!
-					 * mlfan_curr_index is the index of mlfan_next here, and mlfan_next is not the real next one
-					 * (i.e. not the future mlfan_curr)...
-					 */
-					mlfan_curr_index = (e2lfan_curr[0] == mlfan_curr_index) ? e2lfan_curr[1] : e2lfan_curr[0];
-					mpfan_curr_index = loop_to_poly[mlfan_curr_index];
-
-					BLI_assert(mlfan_curr_index >= 0);
-					BLI_assert(mpfan_curr_index >= 0);
-
-					mlfan_next = &mloops[mlfan_curr_index];
-					mpfan_next = &mpolys[mpfan_curr_index];
-					if ((mlfan_curr->v == mlfan_next->v && mlfan_curr->v == mv_pivot_index) ||
-					    (mlfan_curr->v != mlfan_next->v && mlfan_curr->v != mv_pivot_index))
-					{
-						/* We need the previous loop, but current one is our vertex's loop. */
-						mlfan_vert_index = mlfan_curr_index;
-						if (--mlfan_curr_index < mpfan_next->loopstart) {
-							mlfan_curr_index = mpfan_next->loopstart + mpfan_next->totloop - 1;
-						}
+#undef INDEX_UNSET
+#undef INDEX_INVALID
+#undef IS_EDGE_SHARP
+
+/**
+ * Compute internal representation of given custom normals (as an array of float[2]).
+ * It also makes sure the mesh matches those custom normals, by setting sharp edges flag as needed to get a
+ * same custom lnor for all loops sharing a same smooth fan.
+ * If use_vertices if true, custom_loopnors is assumed to be per-vertex, not per-loop
+ * (this allows to set whole vert's normals at once, useful in some cases).
+ */
+static void mesh_normals_loop_custom_set(
+        MVert *mverts, const int numVerts, MEdge *medges, const int numEdges,
+        MLoop *mloops, float (*custom_loopnors)[3], const int numLoops,
+        MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        short (*r_clnors_data)[2], const bool use_vertices)
+{
+	/* We *may* make that poor BKE_mesh_normals_loop_split() even more complex by making it handling that
+	 * feature too, would probably be more efficient in absolute.
+	 * However, this function *is not* performance-critical, since it is mostly expected to be called
+	 * by io addons when importing custom normals, and modifier (and perhaps from some editing tools later?).
+	 * So better to keep some simplicity here, and just call BKE_mesh_normals_loop_split() twice!
+	 */
+	MLoopNorSpaceArray lnors_spacearr = {NULL};
+	BLI_bitmap *done_loops = BLI_BITMAP_NEW((size_t)numLoops, __func__);
+	float (*lnors)[3] = MEM_callocN(sizeof(*lnors) * (size_t)numLoops, __func__);
+	int *loop_to_poly = MEM_mallocN(sizeof(int) * (size_t)numLoops, __func__);
+	/* In this case we always consider split nors as ON, and do not want to use angle to define smooth fans! */
+	const bool use_split_normals = true;
+	const float split_angle = (float)M_PI;
+	int i;
+
+	BLI_SMALLSTACK_DECLARE(clnors_data, short *);
+
+	/* Compute current lnor spacearr. */
+	BKE_mesh_normals_loop_split(mverts, numVerts, medges, numEdges, mloops, lnors, numLoops,
+	                            mpolys, polynors, numPolys, use_split_normals, split_angle,
+	                            &lnors_spacearr, NULL, loop_to_poly);
+
+	/* Now, check each current smooth fan (one lnor space per smooth fan!), and if all its matching custom lnors
+	 * are not (enough) equal, add sharp edges as needed.
+	 * This way, next time we run BKE_mesh_normals_loop_split(), we'll get lnor spacearr/smooth fans matching
+	 * given custom lnors.
+	 * Note this code *will never* unsharp edges!
+	 * And quite obviously, when we set custom normals per vertices, running this is absolutely useless.
+	 */
+	if (!use_vertices) {
+		for (i = 0; i < numLoops; i++) {
+			if (!lnors_spacearr.lspacearr[i]) {
+				/* This should not happen in theory, but in some rare case (probably ugly geometry)
+				 * we can get some NULL loopspacearr at this point. :/
+				 * Maybe we should set those loops' edges as sharp?
+				 */
+				BLI_BITMAP_ENABLE(done_loops, i);
+				printf("WARNING! Getting invalid NULL loop space for loop %d!\n", i);
+				continue;
+			}
+
+			if (!BLI_BITMAP_TEST(done_loops, i)) {
+				/* Notes:
+				 *     * In case of mono-loop smooth fan, loops is NULL, so everything is fine (we have nothing to do).
+				 *     * Loops in this linklist are ordered (in reversed order compared to how they were discovered by
+				 *       BKE_mesh_normals_loop_split(), but this is not a problem). Which means if we find a
+				 *       mismatching clnor, we know all remaining loops will have to be in a new, different smooth fan/
+				 *       lnor space.
+				 *     * In smooth fan case, we compare each clnor against a ref one, to avoid small differences adding
+				 *       up into a real big one in the end!
+				 */
+				LinkNode *loops = lnors_spacearr.lspacearr[i]->loops;
+				MLoop *prev_ml = NULL;
+				const float *org_nor = NULL;
+
+				while (loops) {
+					const int lidx = GET_INT_FROM_POINTER(loops->link);
+					MLoop *ml = &mloops[lidx];
+					const int nidx = use_vertices ? (int)ml->v : lidx;
+					float *nor = custom_loopnors[nidx];
+
+					if (!org_nor) {
+						org_nor = nor;
 					}
-					else {
-						/* We need the next loop, which is also our vertex's loop. */
-						if (++mlfan_curr_index >= mpfan_next->loopstart + mpfan_next->totloop) {
-							mlfan_curr_index = mpfan_next->loopstart;
-						}
-						mlfan_vert_index = mlfan_curr_index;
+					else if (dot_v3v3(org_nor, nor) < LNOR_SPACE_TRIGO_THRESHOLD) {
+						/* Current normal differs too much from org one, we have to tag the edge between
+						 * previous loop's face and current's one as sharp.
+						 * We know those two loops do not point to the same edge, since we do not allow reversed winding
+						 * in a same smooth fan.
+						 */
+						const MPoly *mp = &mpolys[loop_to_poly[lidx]];
+						const MLoop *mlp = &mloops[(lidx == mp->loopstart) ? mp->loopstart + mp->totloop - 1 : lidx - 1];
+						medges[(prev_ml->e == mlp->e) ? prev_ml->e : ml->e].flag |= ME_SHARP;
+
+						org_nor = nor;
 					}
-					mlfan_curr = &mloops[mlfan_curr_index];
-					/* And now we are back in sync, mlfan_curr_index is the index of mlfan_curr! Pff! */
 
-					e2lfan_curr = edge_to_loops[mlfan_curr->e];
+					prev_ml = ml;
+					loops = loops->next;
+					BLI_BITMAP_ENABLE(done_loops, lidx);
 				}
+				BLI_BITMAP_ENABLE(done_loops, i);  /* For single loops, where lnors_spacearr.lspacearr[i]->loops is NULL. */
+			}
+		}
 
-				/* In case we get a zero normal here, just use vertex normal already set! */
-				if (LIKELY(normalize_v3(lnor) != 0.0f)) {
-					/* Copy back the final computed normal into all related loop-normals. */
-					float *nor;
-					while ((nor = BLI_SMALLSTACK_POP(normal))) {
-						copy_v3_v3(nor, lnor);
-					}
+		/* And now, recompute our new auto lnors and lnor spacearr! */
+		BKE_lnor_spacearr_clear(&lnors_spacearr);
+		BKE_mesh_normals_loop_split(mverts, numVerts, medges, numEdges, mloops, lnors, numLoops,
+		                            mpolys, polynors, numPolys, use_split_normals, split_angle,
+		                            &lnors_spacearr, NULL, loop_to_poly);
+	}
+	else {
+		BLI_BITMAP_SET_ALL(done_loops, true, (size_t)numLoops);
+	}
+
+	/* And we just have to convert plain object-space custom normals to our lnor space-encoded ones. */
+	for (i = 0; i < numLoops; i++) {
+		if (!lnors_spacearr.lspacearr[i]) {
+			BLI_BITMAP_DISABLE(done_loops, i);
+			printf("WARNING! Still getting invalid NULL loop space in second loop for loop %d!\n", i);
+			continue;
+		}
+
+		if (BLI_BITMAP_TEST_BOOL(done_loops, i)) {
+			/* Note we accumulate and average all custom normals in current smooth fan, to avoid getting different
+			 * clnors data (tiny differences in plain custom normals can give rather huge differences in
+			 * computed 2D factors).
+			 */
+			LinkNode *loops = lnors_spacearr.lspacearr[i]->loops;
+			if (loops) {
+				int nbr_nors = 0;
+				float avg_nor[3];
+				short clnor_data_tmp[2], *clnor_data;
+
+				zero_v3(avg_nor);
+				while (loops) {
+					const int lidx = GET_INT_FROM_POINTER(loops->link);
+					const int nidx = use_vertices ? (int)mloops[lidx].v : lidx;
+					float *nor = custom_loopnors[nidx];
+
+					nbr_nors++;
+					add_v3_v3(avg_nor, nor);
+					BLI_SMALLSTACK_PUSH(clnors_data, (short *)r_clnors_data[lidx]);
+
+					loops = loops->next;
+					BLI_BITMAP_DISABLE(done_loops, lidx);
 				}
-				else {
-					/* We still have to clear the stack! */
-					while (BLI_SMALLSTACK_POP(normal));
+
+				mul_v3_fl(avg_nor, 1.0f / (float)nbr_nors);
+				BKE_lnor_space_custom_normal_to_data(lnors_spacearr.lspacearr[i], avg_nor, clnor_data_tmp);
+
+				while ((clnor_data = BLI_SMALLSTACK_POP(clnors_data))) {
+					clnor_data[0] = clnor_data_tmp[0];
+					clnor_data[1] = clnor_data_tmp[1];
 				}
 			}
+			else {
+				const int nidx = use_vertices ? (int)mloops[i].v : i;
+				float *nor = custom_loopnors[nidx];
 
-			ml_prev = ml_curr;
-			ml_prev_index = ml_curr_index;
+				BKE_lnor_space_custom_normal_to_data(lnors_spacearr.lspacearr[i], nor, r_clnors_data[i]);
+				BLI_BITMAP_DISABLE(done_loops, i);
+			}
 		}
 	}
 
-	MEM_freeN(edge_to_loops);
+	MEM_freeN(lnors);
 	MEM_freeN(loop_to_poly);
+	MEM_freeN(done_loops);
+	BKE_lnor_spacearr_free(&lnors_spacearr);
+}
 
-#ifdef DEBUG_TIME
-	TIMEIT_END(BKE_mesh_normals_loop_split);
-#endif
-
-#undef INDEX_UNSET
-#undef INDEX_INVALID
-#undef IS_EDGE_SHARP
+void BKE_mesh_normals_loop_custom_set(
+        MVert *mverts, const int numVerts, MEdge *medges, const int numEdges,
+        MLoop *mloops, float (*custom_loopnors)[3], const int numLoops,
+        MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        short (*r_clnors_data)[2])
+{
+	mesh_normals_loop_custom_set(mverts, numVerts, medges, numEdges, mloops, custom_loopnors, numLoops,
+	                             mpolys, polynors, numPolys, r_clnors_data, false);
+}
 
-	}
+void BKE_mesh_normals_loop_custom_from_vertices_set(
+        MVert *mverts, float (*custom_vertnors)[3], const int numVerts,
+        MEdge *medges, const int numEdges, MLoop *mloops, const int numLoops,
+        MPoly *mpolys, const float (*polynors)[3], const int numPolys,
+        short (*r_clnors_data)[2])
+{
+	mesh_normals_loop_custom_set(mverts, numVerts, medges, numEdges, mloops, custom_vertnors, numLoops,
+	                             mpolys, polynors, numPolys, r_clnors_data, true);
 }
 
+#undef LNOR_SPACE_TRIGO_THRESHOLD
 
 /** \} */
 
diff --git a/source/blender/blenkernel/intern/mesh_remap.c b/source/blender/blenkernel/intern/mesh_remap.c
index 7e94e17d1ff..aca72614094 100644
--- a/source/blender/blenkernel/intern/mesh_remap.c
+++ b/source/blender/blenkernel/intern/mesh_remap.c
@@ -1047,6 +1047,8 @@ void BKE_mesh_remap_calc_loops_from_dm(
 				}
 			}
 			if (need_lnors_dst) {
+				short (*custom_nors_dst)[2] = CustomData_get_layer(ldata_dst, CD_CUSTOMLOOPNORMAL);
+
 				/* Cache poly nors into a temp CDLayer. */
 				loop_nors_dst = CustomData_get_layer(ldata_dst, CD_NORMAL);
 				if (dirty_nors_dst || !loop_nors_dst) {
@@ -1056,8 +1058,8 @@ void BKE_mesh_remap_calc_loops_from_dm(
 					}
 					BKE_mesh_normals_loop_split(verts_dst, numverts_dst, edges_dst, numedges_dst,
 					                            loops_dst, loop_nors_dst, numloops_dst,
-					                            polys_dst, poly_nors_dst, numpolys_dst,
-					                            use_split_nors_dst, split_angle_dst);
+					                            polys_dst, (const float (*)[3])poly_nors_dst, numpolys_dst,
+					                            use_split_nors_dst, split_angle_dst, NULL, custom_nors_dst, NULL);
 				}
 			}
 			if (need_pnors_src || need_lnors_src) {
diff --git a/source/blender/blenkernel/intern/subsurf_ccg.c b/source/blender/blenkernel/intern/subsurf_ccg.c
index 45ec3370672..9519c7b25a1 100644
--- a/source/blender/blenkernel/intern/subsurf_ccg.c
+++ b/source/blender/blenkernel/intern/subsurf_ccg.c
@@ -3465,6 +3465,7 @@ static CCGDerivedMesh *getCCGDerivedMesh(CCGSubSurf *ss,
 
 	ccgdm->dm.calcNormals = ccgDM_calcNormals;
 	ccgdm->dm.calcLoopNormals = CDDM_calc_loop_normals;
+	ccgdm->dm.calcLoopNormalsSpaceArray = CDDM_calc_loop_normals_spacearr;
 	ccgdm->dm.recalcTessellation = ccgDM_recalcTessellation;
 
 	ccgdm->dm.getVertCos = ccgdm_getVertCos;
diff --git a/source/blender/blenloader/intern/versioning_270.c b/source/blender/blenloader/intern/versioning_270.c
index 096fc315f39..4ca1a44d64e 100644
--- a/source/blender/blenloader/intern/versioning_270.c
+++ b/source/blender/blenloader/intern/versioning_270.c
@@ -473,90 +473,92 @@ void blo_do_versions_270(FileData *fd, Library *UNUSED(lib), Main *main)
 		}
 	}
 	
-	if (!DNA_struct_elem_find(fd->filesdna, "ClothSimSettings", "float", "bending_damping")) {
-		Object *ob;
-		ModifierData *md;
-		for (ob = main->object.first; ob; ob = ob->id.next) {
-			for (md = ob->modifiers.first; md; md = md->next) {
-				if (md->type == eModifierType_Cloth) {
-					ClothModifierData *clmd = (ClothModifierData *)md;
-					clmd->sim_parms->bending_damping = 0.5f;
-				}
-				else if (md->type == eModifierType_ParticleSystem) {
-					ParticleSystemModifierData *pmd = (ParticleSystemModifierData *)md;
-					if (pmd->psys->clmd) {
-						pmd->psys->clmd->sim_parms->bending_damping = 0.5f;
+	if (!MAIN_VERSION_ATLEAST(main, 273, 3)) {
+		ParticleSettings *part;
+		for (part = main->particle.first; part; part = part->id.next) {
+			if (part->clumpcurve)
+				part->child_flag |= PART_CHILD_USE_CLUMP_CURVE;
+			if (part->roughcurve)
+				part->child_flag |= PART_CHILD_USE_ROUGH_CURVE;
+		}
+	}
+
+	if (!MAIN_VERSION_ATLEAST(main, 273, 6)) {
+		if (!DNA_struct_elem_find(fd->filesdna, "ClothSimSettings", "float", "bending_damping")) {
+			Object *ob;
+			ModifierData *md;
+			for (ob = main->object.first; ob; ob = ob->id.next) {
+				for (md = ob->modifiers.first; md; md = md->next) {
+					if (md->type == eModifierType_Cloth) {
+						ClothModifierData *clmd = (ClothModifierData *)md;
+						clmd->sim_parms->bending_damping = 0.5f;
+					}
+					else if (md->type == eModifierType_ParticleSystem) {
+						ParticleSystemModifierData *pmd = (ParticleSystemModifierData *)md;
+						if (pmd->psys->clmd) {
+							pmd->psys->clmd->sim_parms->bending_damping = 0.5f;
+						}
 					}
 				}
 			}
 		}
-	}
 
-	if (!DNA_struct_elem_find(fd->filesdna, "ParticleSettings", "float", "clump_noise_size")) {
-		ParticleSettings *part;
-		for (part = main->particle.first; part; part = part->id.next) {
-			part->clump_noise_size = 1.0f;
+		if (!DNA_struct_elem_find(fd->filesdna, "ParticleSettings", "float", "clump_noise_size")) {
+			ParticleSettings *part;
+			for (part = main->particle.first; part; part = part->id.next) {
+				part->clump_noise_size = 1.0f;
+			}
 		}
-	}
 
-	if (!DNA_struct_elem_find(fd->filesdna, "ParticleSettings", "int", "kink_extra_steps")) {
-		ParticleSettings *part;
-		for (part = main->particle.first; part; part = part->id.next) {
-			part->kink_extra_steps = 4;
+		if (!DNA_struct_elem_find(fd->filesdna, "ParticleSettings", "int", "kink_extra_steps")) {
+			ParticleSettings *part;
+			for (part = main->particle.first; part; part = part->id.next) {
+				part->kink_extra_steps = 4;
+			}
 		}
-	}
 
-	if (!DNA_struct_elem_find(fd->filesdna, "MTex", "float", "kinkampfac")) {
-		ParticleSettings *part;
-		for (part = main->particle.first; part; part = part->id.next) {
-			int a;
-			for (a = 0; a < MAX_MTEX; a++) {
-				MTex *mtex = part->mtex[a];
-				if (mtex) {
-					mtex->kinkampfac = 1.0f;
+		if (!DNA_struct_elem_find(fd->filesdna, "MTex", "float", "kinkampfac")) {
+			ParticleSettings *part;
+			for (part = main->particle.first; part; part = part->id.next) {
+				int a;
+				for (a = 0; a < MAX_MTEX; a++) {
+					MTex *mtex = part->mtex[a];
+					if (mtex) {
+						mtex->kinkampfac = 1.0f;
+					}
 				}
 			}
 		}
-	}
 
-	if (!DNA_struct_elem_find(fd->filesdna, "HookModifierData", "char", "flag")) {
-		Object *ob;
+		if (!DNA_struct_elem_find(fd->filesdna, "HookModifierData", "char", "flag")) {
+			Object *ob;
 
-		for (ob = main->object.first; ob; ob = ob->id.next) {
-			ModifierData *md;
-			for (md = ob->modifiers.first; md; md = md->next) {
-				if (md->type == eModifierType_Hook) {
-					HookModifierData *hmd = (HookModifierData *)md;
-					hmd->falloff_type = eHook_Falloff_InvSquare;
+			for (ob = main->object.first; ob; ob = ob->id.next) {
+				ModifierData *md;
+				for (md = ob->modifiers.first; md; md = md->next) {
+					if (md->type == eModifierType_Hook) {
+						HookModifierData *hmd = (HookModifierData *)md;
+						hmd->falloff_type = eHook_Falloff_InvSquare;
+					}
 				}
 			}
 		}
-	}
 
-	if (!MAIN_VERSION_ATLEAST(main, 273, 3)) {
-		ParticleSettings *part;
-		for (part = main->particle.first; part; part = part->id.next) {
-			if (part->clumpcurve)
-				part->child_flag |= PART_CHILD_USE_CLUMP_CURVE;
-			if (part->roughcurve)
-				part->child_flag |= PART_CHILD_USE_ROUGH_CURVE;
-		}
-	}
-
-	if (!DNA_struct_elem_find(fd->filesdna, "NodePlaneTrackDeformData", "char", "flag")) {
-		FOREACH_NODETREE(main, ntree, id) {
-			if (ntree->type == NTREE_COMPOSIT) {
-				bNode *node;
-				for (node = ntree->nodes.first; node; node = node->next) {
-					if (ELEM(node->type, CMP_NODE_PLANETRACKDEFORM)) {
-						NodePlaneTrackDeformData *data = node->storage;
-						data->flag = 0;
-						data->motion_blur_samples = 16;
-						data->motion_blur_shutter = 0.5f;
+		if (!DNA_struct_elem_find(fd->filesdna, "NodePlaneTrackDeformData", "char", "flag")) {
+			FOREACH_NODETREE(main, ntree, id) {
+				if (ntree->type == NTREE_COMPOSIT) {
+					bNode *node;
+					for (node = ntree->nodes.first; node; node = node->next) {
+						if (ELEM(node->type, CMP_NODE_PLANETRACKDEFORM)) {
+							NodePlaneTrackDeformData *data = node->storage;
+							data->flag = 0;
+							data->motion_blur_samples = 16;
+							data->motion_blur_shutter = 0.5f;
+						}
 					}
 				}
 			}
+			FOREACH_NODETREE_END
 		}
-		FOREACH_NODETREE_END
 	}
 }
diff --git a/source/blender/bmesh/intern/bmesh_mesh.c b/source/blender/bmesh/intern/bmesh_mesh.c
index 95d40bfc0c6..6700ffe314e 100644
--- a/source/blender/bmesh/intern/bmesh_mesh.c
+++ b/source/blender/bmesh/intern/bmesh_mesh.c
@@ -34,10 +34,12 @@
 #include "BLI_linklist_stack.h"
 #include "BLI_listbase.h"
 #include "BLI_math.h"
+#include "BLI_stack.h"
 #include "BLI_utildefines.h"
 
 #include "BKE_cdderivedmesh.h"
 #include "BKE_editmesh.h"
+#include "BKE_mesh.h"
 #include "BKE_multires.h"
 
 #include "intern/bmesh_private.h"
@@ -438,7 +440,8 @@ void BM_verts_calc_normal_vcos(BMesh *bm, const float (*fnos)[3], const float (*
 static void bm_mesh_edges_sharp_tag(BMesh *bm, const float (*vnos)[3], const float (*fnos)[3], float split_angle,
                                     float (*r_lnos)[3])
 {
-	BMIter eiter;
+	BMIter eiter, viter;
+	BMVert *v;
 	BMEdge *e;
 	int i;
 
@@ -450,15 +453,18 @@ static void bm_mesh_edges_sharp_tag(BMesh *bm, const float (*vnos)[3], const flo
 
 	{
 		char htype = BM_LOOP;
-		if (vnos) {
-			htype |= BM_VERT;
-		}
 		if (fnos) {
 			htype |= BM_FACE;
 		}
 		BM_mesh_elem_index_ensure(bm, htype);
 	}
 
+	/* Clear all vertices' tags (means they are all smooth for now). */
+	BM_ITER_MESH_INDEX (v, &viter, bm, BM_VERTS_OF_MESH, i) {
+		BM_elem_index_set(v, i); /* set_inline */
+		BM_elem_flag_disable(v, BM_ELEM_TAG);
+	}
+
 	/* This first loop checks which edges are actually smooth, and pre-populate lnos with vnos (as if they were
 	 * all smooth).
 	 */
@@ -481,13 +487,13 @@ static void bm_mesh_edges_sharp_tag(BMesh *bm, const float (*vnos)[3], const flo
 			 * If the angle between both its polys' normals is below split_angle value,
 			 * and it is tagged as such,
 			 * and both its faces are smooth,
-			 * and both its faces have compatible (non-flipped) normals, i.e. both loops on the same edge do not share
-			 *     the same vertex.
+			 * and both its faces have compatible (non-flipped) normals,
+			 * i.e. both loops on the same edge do not share the same vertex.
 			 */
 			if (is_angle_smooth &&
-			    BM_elem_flag_test_bool(e, BM_ELEM_SMOOTH) &&
-			    BM_elem_flag_test_bool(l_a->f, BM_ELEM_SMOOTH) &&
-			    BM_elem_flag_test_bool(l_b->f, BM_ELEM_SMOOTH) &&
+			    BM_elem_flag_test(e, BM_ELEM_SMOOTH) &&
+			    BM_elem_flag_test(l_a->f, BM_ELEM_SMOOTH) &&
+			    BM_elem_flag_test(l_b->f, BM_ELEM_SMOOTH) &&
 			    l_a->v != l_b->v)
 			{
 				const float *no;
@@ -499,20 +505,40 @@ static void bm_mesh_edges_sharp_tag(BMesh *bm, const float (*vnos)[3], const flo
 				no = vnos ? vnos[BM_elem_index_get(l_b->v)] : l_b->v->no;
 				copy_v3_v3(r_lnos[BM_elem_index_get(l_b)], no);
 			}
+			else {
+				/* Sharp edge, tag its verts as such. */
+				BM_elem_flag_enable(e->v1, BM_ELEM_TAG);
+				BM_elem_flag_enable(e->v2, BM_ELEM_TAG);
+			}
+		}
+		else {
+			/* Sharp edge, tag its verts as such. */
+			BM_elem_flag_enable(e->v1, BM_ELEM_TAG);
+			BM_elem_flag_enable(e->v2, BM_ELEM_TAG);
 		}
 	}
 
-	bm->elem_index_dirty &= ~BM_EDGE;
+	bm->elem_index_dirty &= ~(BM_EDGE | BM_VERT);
 }
 
-/* BMesh version of BKE_mesh_normals_loop_split() in mesh_evaluate.c */
-static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const float (*fnos)[3], float (*r_lnos)[3])
+/* BMesh version of BKE_mesh_normals_loop_split() in mesh_evaluate.c
+ * Will use first clnors_data array, and fallback to cd_loop_clnors_offset (use NULL and -1 to not use clnors). */
+static void bm_mesh_loops_calc_normals(
+        BMesh *bm, const float (*vcos)[3], const float (*fnos)[3], float (*r_lnos)[3],
+        MLoopNorSpaceArray *r_lnors_spacearr, short (*clnors_data)[2], const int cd_loop_clnors_offset)
 {
 	BMIter fiter;
 	BMFace *f_curr;
+	const bool has_clnors = clnors_data || (cd_loop_clnors_offset != -1);
+
+	MLoopNorSpaceArray _lnors_spacearr = {NULL};
 
 	/* Temp normal stack. */
 	BLI_SMALLSTACK_DECLARE(normal, float *);
+	/* Temp clnors stack. */
+	BLI_SMALLSTACK_DECLARE(clnors, short *);
+	/* Temp edge vectors stack, only used when computing lnor spacearr. */
+	BLI_Stack *edge_vectors = NULL;
 
 	{
 		char htype = BM_LOOP;
@@ -525,6 +551,15 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 		BM_mesh_elem_index_ensure(bm, htype);
 	}
 
+	if (!r_lnors_spacearr && has_clnors) {
+		/* We need to compute lnor spacearr if some custom lnor data are given to us! */
+		r_lnors_spacearr = &_lnors_spacearr;
+	}
+	if (r_lnors_spacearr) {
+		BKE_lnor_spacearr_init(r_lnors_spacearr, bm->totloop);
+		edge_vectors = BLI_stack_new(sizeof(float[3]), __func__);
+	}
+
 	/* We now know edges that can be smoothed (they are tagged), and edges that will be hard (they aren't).
 	 * Now, time to generate the normals.
 	 */
@@ -533,8 +568,10 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 
 		l_curr = l_first = BM_FACE_FIRST_LOOP(f_curr);
 		do {
-			if (BM_elem_flag_test_bool(l_curr->e, BM_ELEM_TAG)) {
-				/* A smooth edge.
+			if (BM_elem_flag_test(l_curr->e, BM_ELEM_TAG) &&
+			    (!r_lnors_spacearr || BM_elem_flag_test(l_curr->v, BM_ELEM_TAG)))
+			{
+				/* A smooth edge, and we are not generating lnors_spacearr, or the related vertex is sharp.
 				 * We skip it because it is either:
 				 * - in the middle of a 'smooth fan' already computed (or that will be as soon as we hit
 				 *   one of its ends, i.e. one of its two sharp edges), or...
@@ -542,12 +579,45 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 				 *   are just fine!
 				 */
 			}
-			else if (!BM_elem_flag_test_bool(l_curr->prev->e, BM_ELEM_TAG)) {
+			else if (!BM_elem_flag_test(l_curr->e, BM_ELEM_TAG) &&
+			         !BM_elem_flag_test(l_curr->prev->e, BM_ELEM_TAG))
+			{
 				/* Simple case (both edges around that vertex are sharp in related polygon),
 				 * this vertex just takes its poly normal.
 				 */
+				const int l_curr_index = BM_elem_index_get(l_curr);
 				const float *no = fnos ? fnos[BM_elem_index_get(f_curr)] : f_curr->no;
-				copy_v3_v3(r_lnos[BM_elem_index_get(l_curr)], no);
+				copy_v3_v3(r_lnos[l_curr_index], no);
+
+				/* If needed, generate this (simple!) lnor space. */
+				if (r_lnors_spacearr) {
+					float vec_curr[3], vec_prev[3];
+					MLoopNorSpace *lnor_space = BKE_lnor_space_create(r_lnors_spacearr);
+
+					{
+						const BMVert *v_pivot = l_curr->v;
+						const float *co_pivot = vcos ? vcos[BM_elem_index_get(v_pivot)] : v_pivot->co;
+						const BMVert *v_1 = BM_edge_other_vert(l_curr->e, v_pivot);
+						const float *co_1 = vcos ? vcos[BM_elem_index_get(v_1)] : v_1->co;
+						const BMVert *v_2 = BM_edge_other_vert(l_curr->prev->e, v_pivot);
+						const float *co_2 = vcos ? vcos[BM_elem_index_get(v_2)] : v_2->co;
+
+						sub_v3_v3v3(vec_curr, co_1, co_pivot);
+						normalize_v3(vec_curr);
+						sub_v3_v3v3(vec_prev, co_2, co_pivot);
+						normalize_v3(vec_prev);
+					}
+
+					BKE_lnor_space_define(lnor_space, r_lnos[l_curr_index], vec_curr, vec_prev, NULL);
+					/* We know there is only one loop in this space, no need to create a linklist in this case... */
+					BKE_lnor_space_add_loop(r_lnors_spacearr, lnor_space, l_curr_index, false);
+
+					if (has_clnors) {
+						short (*clnor)[2] = clnors_data ? &clnors_data[l_curr_index] :
+						                                  BM_ELEM_CD_GET_VOID_P(l_curr, cd_loop_clnors_offset);
+						BKE_lnor_space_custom_data_to_normal(lnor_space, *clnor, r_lnos[l_curr_index]);
+					}
+				}
 			}
 			/* We *do not need* to check/tag loops as already computed!
 			 * Due to the fact a loop only links to one of its two edges, a same fan *will never be walked more than
@@ -567,13 +637,26 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 				 */
 				BMVert *v_pivot = l_curr->v;
 				BMEdge *e_next;
+				const BMEdge *e_org = l_curr->e;
 				BMLoop *lfan_pivot, *lfan_pivot_next;
+				int lfan_pivot_index;
 				float lnor[3] = {0.0f, 0.0f, 0.0f};
-				float vec_curr[3], vec_next[3];
+				float vec_curr[3], vec_next[3], vec_org[3];
+
+				/* We validate clnors data on the fly - cheapest way to do! */
+				int clnors_avg[2] = {0, 0};
+				short (*clnor_ref)[2] = NULL;
+				int clnors_nbr = 0;
+				bool clnors_invalid = false;
 
 				const float *co_pivot = vcos ? vcos[BM_elem_index_get(v_pivot)] : v_pivot->co;
 
+				MLoopNorSpace *lnor_space = r_lnors_spacearr ? BKE_lnor_space_create(r_lnors_spacearr) : NULL;
+
+				BLI_assert((edge_vectors == NULL) || BLI_stack_is_empty(edge_vectors));
+
 				lfan_pivot = l_curr;
+				lfan_pivot_index = BM_elem_index_get(lfan_pivot);
 				e_next = lfan_pivot->e;  /* Current edge here, actually! */
 
 				/* Only need to compute previous edge's vector once, then we can just reuse old current one! */
@@ -581,8 +664,13 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 					const BMVert *v_2 = BM_edge_other_vert(e_next, v_pivot);
 					const float *co_2 = vcos ? vcos[BM_elem_index_get(v_2)] : v_2->co;
 
-					sub_v3_v3v3(vec_curr, co_2, co_pivot);
-					normalize_v3(vec_curr);
+					sub_v3_v3v3(vec_org, co_2, co_pivot);
+					normalize_v3(vec_org);
+					copy_v3_v3(vec_curr, vec_org);
+
+					if (r_lnors_spacearr) {
+						BLI_stack_push(edge_vectors, vec_org);
+					}
 				}
 
 				while (true) {
@@ -617,12 +705,38 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 						const float *no = fnos ? fnos[BM_elem_index_get(f)] : f->no;
 						/* Accumulate */
 						madd_v3_v3fl(lnor, no, fac);
+
+						if (has_clnors) {
+							/* Accumulate all clnors, if they are not all equal we have to fix that! */
+							short (*clnor)[2] = clnors_data ? &clnors_data[lfan_pivot_index] :
+							                                  BM_ELEM_CD_GET_VOID_P(lfan_pivot, cd_loop_clnors_offset);
+							if (clnors_nbr) {
+								clnors_invalid |= ((*clnor_ref)[0] != (*clnor)[0] || (*clnor_ref)[1] != (*clnor)[1]);
+							}
+							else {
+								clnor_ref = clnor;
+							}
+							clnors_avg[0] += (*clnor)[0];
+							clnors_avg[1] += (*clnor)[1];
+							clnors_nbr++;
+							/* We store here a pointer to all custom lnors processed. */
+							BLI_SMALLSTACK_PUSH(clnors, (short *)*clnor);
+						}
 					}
 
 					/* We store here a pointer to all loop-normals processed. */
-					BLI_SMALLSTACK_PUSH(normal, (float *)r_lnos[BM_elem_index_get(lfan_pivot)]);
+					BLI_SMALLSTACK_PUSH(normal, (float *)r_lnos[lfan_pivot_index]);
+
+					if (r_lnors_spacearr) {
+						/* Assign current lnor space to current 'vertex' loop. */
+						BKE_lnor_space_add_loop(r_lnors_spacearr, lnor_space, lfan_pivot_index, true);
+						if (e_next != e_org) {
+							/* We store here all edges-normalized vectors processed. */
+							BLI_stack_push(edge_vectors, vec_next);
+						}
+					}
 
-					if (!BM_elem_flag_test_bool(e_next, BM_ELEM_TAG)) {
+					if (!BM_elem_flag_test(e_next, BM_ELEM_TAG) || (e_next == e_org)) {
 						/* Next edge is sharp, we have finished with this fan of faces around this vert! */
 						break;
 					}
@@ -631,23 +745,74 @@ static void bm_mesh_loops_calc_normals(BMesh *bm, const float (*vcos)[3], const
 					copy_v3_v3(vec_curr, vec_next);
 					/* Next pivot loop to current one. */
 					lfan_pivot = lfan_pivot_next;
+					lfan_pivot_index = BM_elem_index_get(lfan_pivot);
 				}
 
-				/* In case we get a zero normal here, just use vertex normal already set! */
-				if (LIKELY(normalize_v3(lnor) != 0.0f)) {
-					/* Copy back the final computed normal into all related loop-normals. */
-					float *nor;
-					while ((nor = BLI_SMALLSTACK_POP(normal))) {
-						copy_v3_v3(nor, lnor);
+				{
+					float lnor_len = normalize_v3(lnor);
+
+					/* If we are generating lnor spacearr, we can now define the one for this fan. */
+					if (r_lnors_spacearr) {
+						if (UNLIKELY(lnor_len == 0.0f)) {
+							/* Use vertex normal as fallback! */
+							copy_v3_v3(lnor, r_lnos[lfan_pivot_index]);
+							lnor_len = 1.0f;
+						}
+
+						BKE_lnor_space_define(lnor_space, lnor, vec_org, vec_next, edge_vectors);
+
+						if (has_clnors) {
+							if (clnors_invalid) {
+								short *clnor;
+
+								clnors_avg[0] /= clnors_nbr;
+								clnors_avg[1] /= clnors_nbr;
+								/* Fix/update all clnors of this fan with computed average value. */
+								printf("Invalid clnors in this fan!\n");
+								while ((clnor = BLI_SMALLSTACK_POP(clnors))) {
+									//print_v2("org clnor", clnor);
+									clnor[0] = (short)clnors_avg[0];
+									clnor[1] = (short)clnors_avg[1];
+								}
+								//print_v2("new clnors", clnors_avg);
+							}
+							else {
+								/* We still have to consume the stack! */
+								while (BLI_SMALLSTACK_POP(clnors));
+							}
+							BKE_lnor_space_custom_data_to_normal(lnor_space, *clnor_ref, lnor);
+						}
+					}
+
+					/* In case we get a zero normal here, just use vertex normal already set! */
+					if (LIKELY(lnor_len != 0.0f)) {
+						/* Copy back the final computed normal into all related loop-normals. */
+						float *nor;
+
+						while ((nor = BLI_SMALLSTACK_POP(normal))) {
+							copy_v3_v3(nor, lnor);
+						}
+					}
+					else {
+						/* We still have to consume the stack! */
+						while (BLI_SMALLSTACK_POP(normal));
 					}
 				}
-				else {
-					/* We still have to clear the stack! */
-					while (BLI_SMALLSTACK_POP(normal));
+
+				/* Tag related vertex as sharp, to avoid fanning around it again (in case it was a smooth one). */
+				if (r_lnors_spacearr) {
+					BM_elem_flag_enable(l_curr->v, BM_ELEM_TAG);
 				}
 			}
 		} while ((l_curr = l_curr->next) != l_first);
 	}
+
+	if (r_lnors_spacearr) {
+		BLI_stack_free(edge_vectors);
+		if (r_lnors_spacearr == &_lnors_spacearr) {
+			BKE_lnor_spacearr_free(r_lnors_spacearr);
+		}
+	}
 }
 
 static void bm_mesh_loops_from_vert_normals(BMesh *bm, const float (*vnos)[3], float (*r_lnos)[3])
@@ -682,16 +847,22 @@ static void bm_mesh_loops_from_vert_normals(BMesh *bm, const float (*vnos)[3], f
  * Updates the loop normals of a mesh. Assumes vertex and face normals are valid (else call BM_mesh_normals_update()
  * first)!
  */
-void BM_mesh_loop_normals_update(BMesh *bm, const bool use_split_normals, const float split_angle, float (*r_lnos)[3])
+void BM_mesh_loop_normals_update(
+        BMesh *bm, const bool use_split_normals, const float split_angle, float (*r_lnos)[3],
+        MLoopNorSpaceArray *r_lnors_spacearr, short (*clnors_data)[2], const int cd_loop_clnors_offset)
 {
+	const bool has_clnors = clnors_data || (cd_loop_clnors_offset != -1);
+
 	if (use_split_normals) {
-		/* Tag smooth edges and set lnos from vnos when they might be completely smooth... */
-		bm_mesh_edges_sharp_tag(bm, NULL, NULL, split_angle, r_lnos);
+	        /* Tag smooth edges and set lnos from vnos when they might be completely smooth...
+	         * When using custom loop normals, disable the angle feature! */
+	        bm_mesh_edges_sharp_tag(bm, NULL, NULL, has_clnors ? (float)M_PI : split_angle, r_lnos);
 
-		/* Finish computing lnos by accumulating face normals in each fan of faces defined by sharp edges. */
-		bm_mesh_loops_calc_normals(bm, NULL, NULL, r_lnos);
+	        /* Finish computing lnos by accumulating face normals in each fan of faces defined by sharp edges. */
+	        bm_mesh_loops_calc_normals(bm, NULL, NULL, r_lnos, r_lnors_spacearr, clnors_data, cd_loop_clnors_offset);
 	}
 	else {
+		BLI_assert(!r_lnors_spacearr);
 		bm_mesh_loops_from_vert_normals(bm, NULL, r_lnos);
 	}
 }
@@ -703,17 +874,23 @@ void BM_mesh_loop_normals_update(BMesh *bm, const bool use_split_normals, const
  * Compute split normals, i.e. vertex normals associated with each poly (hence 'loop normals').
  * Useful to materialize sharp edges (or non-smooth faces) without actually modifying the geometry (splitting edges).
  */
-void BM_loops_calc_normal_vcos(BMesh *bm, const float (*vcos)[3], const float (*vnos)[3], const float (*fnos)[3],
-                               const bool use_split_normals, const float split_angle, float (*r_lnos)[3])
+void BM_loops_calc_normal_vcos(
+        BMesh *bm, const float (*vcos)[3], const float (*vnos)[3], const float (*fnos)[3],
+        const bool use_split_normals, const float split_angle, float (*r_lnos)[3],
+        MLoopNorSpaceArray *r_lnors_spacearr, short (*clnors_data)[2], const int cd_loop_clnors_offset)
 {
+	const bool has_clnors = clnors_data || (cd_loop_clnors_offset != -1);
+
 	if (use_split_normals) {
-		/* Tag smooth edges and set lnos from vnos when they might be completely smooth... */
-		bm_mesh_edges_sharp_tag(bm, vnos, fnos, split_angle, r_lnos);
+		/* Tag smooth edges and set lnos from vnos when they might be completely smooth...
+		 * When using custom loop normals, disable the angle feature! */
+		bm_mesh_edges_sharp_tag(bm, vnos, fnos, has_clnors ? (float)M_PI : split_angle, r_lnos);
 
 		/* Finish computing lnos by accumulating face normals in each fan of faces defined by sharp edges. */
-		bm_mesh_loops_calc_normals(bm, vcos, fnos, r_lnos);
+		bm_mesh_loops_calc_normals(bm, vcos, fnos, r_lnos, r_lnors_spacearr, clnors_data, cd_loop_clnors_offset);
 	}
 	else {
+		BLI_assert(!r_lnors_spacearr);
 		bm_mesh_loops_from_vert_normals(bm, vnos, r_lnos);
 	}
 }
diff --git a/source/blender/bmesh/intern/bmesh_mesh.h b/source/blender/bmesh/intern/bmesh_mesh.h
index c2ac2c40534..bac5da8347e 100644
--- a/source/blender/bmesh/intern/bmesh_mesh.h
+++ b/source/blender/bmesh/intern/bmesh_mesh.h
@@ -28,6 +28,7 @@
  */
 
 struct BMAllocTemplate;
+struct MLoopNorSpaceArray;
 
 void   BM_mesh_elem_toolflags_ensure(BMesh *bm);
 void   BM_mesh_elem_toolflags_clear(BMesh *bm);
@@ -39,8 +40,10 @@ void   BM_mesh_clear(BMesh *bm);
 
 void BM_mesh_normals_update(BMesh *bm);
 void BM_verts_calc_normal_vcos(BMesh *bm, const float (*fnos)[3], const float (*vcos)[3], float (*vnos)[3]);
-void BM_loops_calc_normal_vcos(BMesh *bm, const float (*vcos)[3], const float (*vnos)[3], const float (*pnos)[3],
-                               const bool use_split_normals, const float split_angle, float (*r_lnos)[3]);
+void BM_loops_calc_normal_vcos(
+        BMesh *bm, const float (*vcos)[3], const float (*vnos)[3], const float (*pnos)[3],
+        const bool use_split_normals, const float split_angle, float (*r_lnos)[3],
+        struct MLoopNorSpaceArray *r_lnors_spacearr, short (*clnors_data)[2], const int cd_loop_clnors_offset);
 
 void bmesh_edit_begin(BMesh *bm, const BMOpTypeFlag type_flag);
 void bmesh_edit_end(BMesh *bm, const BMOpTypeFlag type_flag);
diff --git a/source/blender/editors/mesh/mesh_data.c b/source/blender/editors/mesh/mesh_data.c
index 152d055d239..30e9a85d083 100644
--- a/source/blender/editors/mesh/mesh_data.c
+++ b/source/blender/editors/mesh/mesh_data.c
@@ -57,6 +57,7 @@
 
 #include "ED_mesh.h"
 #include "ED_object.h"
+#include "ED_screen.h"
 #include "ED_uvedit.h"
 #include "ED_view3d.h"
 
@@ -838,6 +839,71 @@ void MESH_OT_customdata_clear_skin(wmOperatorType *ot)
 	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
 }
 
+/* Clear custom loop normals */
+static int mesh_customdata_custom_splitnormals_add_exec(bContext *C, wmOperator *UNUSED(op))
+{
+	Object *ob = ED_object_context(C);
+	Mesh *me = ob->data;
+
+	if (!BKE_mesh_has_custom_loop_normals(me)) {
+		CustomData *data = GET_CD_DATA(me, ldata);
+
+		if (me->edit_btmesh) {
+			BM_data_layer_add(me->edit_btmesh->bm, data, CD_CUSTOMLOOPNORMAL);
+		}
+		else {
+			CustomData_add_layer(data, CD_CUSTOMLOOPNORMAL, CD_DEFAULT, NULL, me->totloop);
+		}
+
+		DAG_id_tag_update(&me->id, 0);
+		WM_event_add_notifier(C, NC_GEOM | ND_DATA, me);
+
+		return OPERATOR_FINISHED;
+	}
+	return OPERATOR_CANCELLED;
+}
+
+void MESH_OT_customdata_custom_splitnormals_add(wmOperatorType *ot)
+{
+	/* identifiers */
+	ot->name = "Add Custom Split Normals Data";
+	ot->idname = "MESH_OT_customdata_custom_splitnormals_add";
+	ot->description = "Add a custom split normals layer, if none exists yet";
+
+	/* api callbacks */
+	ot->exec = mesh_customdata_custom_splitnormals_add_exec;
+	ot->poll = ED_operator_object_active_editable_mesh;
+
+	/* flags */
+	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
+}
+
+static int mesh_customdata_custom_splitnormals_clear_exec(bContext *C, wmOperator *UNUSED(op))
+{
+	Object *ob = ED_object_context(C);
+	Mesh *me = ob->data;
+
+	if (BKE_mesh_has_custom_loop_normals(me)) {
+		return mesh_customdata_clear_exec__internal(C, BM_LOOP, CD_CUSTOMLOOPNORMAL);
+	}
+	return OPERATOR_CANCELLED;
+}
+
+void MESH_OT_customdata_custom_splitnormals_clear(wmOperatorType *ot)
+{
+	/* identifiers */
+	ot->name = "Clear Custom Split Normals Data";
+	ot->idname = "MESH_OT_customdata_custom_splitnormals_clear";
+	ot->description = "Remove the custom split normals layer, if it exists";
+
+	/* api callbacks */
+	ot->exec = mesh_customdata_custom_splitnormals_clear_exec;
+	ot->poll = ED_operator_object_active_editable_mesh;
+
+	/* flags */
+	ot->flag = OPTYPE_REGISTER | OPTYPE_UNDO;
+}
+
 /************************** Add Geometry Layers *************************/
 
 void ED_mesh_update(Mesh *mesh, bContext *C, int calc_edges, int calc_tessface)
diff --git a/source/blender/editors/mesh/mesh_intern.h b/source/blender/editors/mesh/mesh_intern.h
index a794db62e3e..8f5ecaed524 100644
--- a/source/blender/editors/mesh/mesh_intern.h
+++ b/source/blender/editors/mesh/mesh_intern.h
@@ -238,6 +238,8 @@ void MESH_OT_vertex_color_remove(struct wmOperatorType *ot);
 /* no create_mask yet */
 void MESH_OT_customdata_clear_mask(struct wmOperatorType *ot);
 void MESH_OT_customdata_clear_skin(struct wmOperatorType *ot);
+void MESH_OT_customdata_custom_splitnormals_add(struct wmOperatorType *ot);
+void MESH_OT_customdata_custom_splitnormals_clear(struct wmOperatorType *ot);
 void MESH_OT_drop_named_image(struct wmOperatorType *ot);
 
 
diff --git a/source/blender/editors/mesh/mesh_ops.c b/source/blender/editors/mesh/mesh_ops.c
index 8146be7ee47..59d0de87de6 100644
--- a/source/blender/editors/mesh/mesh_ops.c
+++ b/source/blender/editors/mesh/mesh_ops.c
@@ -153,6 +153,8 @@ void ED_operatortypes_mesh(void)
 	WM_operatortype_append(MESH_OT_vertex_color_remove);
 	WM_operatortype_append(MESH_OT_customdata_clear_mask);
 	WM_operatortype_append(MESH_OT_customdata_clear_skin);
+	WM_operatortype_append(MESH_OT_customdata_custom_splitnormals_add);
+	WM_operatortype_append(MESH_OT_customdata_custom_splitnormals_clear);
 	WM_operatortype_append(MESH_OT_drop_named_image);
 
 	WM_operatortype_append(MESH_OT_edgering_select);
diff --git a/source/blender/makesdna/DNA_customdata_types.h b/source/blender/makesdna/DNA_customdata_types.h
index 7fedad2546c..74f5967db13 100644
--- a/source/blender/makesdna/DNA_customdata_types.h
+++ b/source/blender/makesdna/DNA_customdata_types.h
@@ -63,9 +63,10 @@ typedef struct CustomDataExternal {
  * layers, each with a data type (e.g. MTFace, MDeformVert, etc.). */
 typedef struct CustomData {
 	CustomDataLayer *layers;      /* CustomDataLayers, ordered by type */
-	int typemap[41];              /* runtime only! - maps types to indices of first layer of that type,
+	int typemap[42];              /* runtime only! - maps types to indices of first layer of that type,
 	                               * MUST be >= CD_NUMTYPES, but we cant use a define here.
 	                               * Correct size is ensured in CustomData_update_typemap assert() */
+	int pad_i1;
 	int totlayer, maxlayer;       /* number of layers, size of layers array */
 	int totsize;                  /* in editmode, total size of all data layers */
 	struct BLI_mempool *pool;     /* (BMesh Only): Memory pool for allocation of blocks */
@@ -119,8 +120,9 @@ typedef enum CustomDataType {
 	CD_FREESTYLE_FACE   = 38,
 	CD_MLOOPTANGENT     = 39,
 	CD_TESSLOOPNORMAL   = 40,
+	CD_CUSTOMLOOPNORMAL = 41,
 
-	CD_NUMTYPES         = 41
+	CD_NUMTYPES         = 42
 } CustomDataType;
 
 /* Bits for CustomDataMask */
@@ -167,6 +169,7 @@ typedef enum CustomDataType {
 #define CD_MASK_FREESTYLE_FACE	(1LL << CD_FREESTYLE_FACE)
 #define CD_MASK_MLOOPTANGENT    (1LL << CD_MLOOPTANGENT)
 #define CD_MASK_TESSLOOPNORMAL  (1LL << CD_TESSLOOPNORMAL)
+#define CD_MASK_CUSTOMLOOPNORMAL (1LL << CD_CUSTOMLOOPNORMAL)
 
 /* CustomData.flag */
 enum {
diff --git a/source/blender/makesrna/intern/rna_mesh.c b/source/blender/makesrna/intern/rna_mesh.c
index 3d6eab2bc88..44bae770186 100644
--- a/source/blender/makesrna/intern/rna_mesh.c
+++ b/source/blender/makesrna/intern/rna_mesh.c
@@ -207,6 +207,11 @@ static void rna_MeshAnyLayer_name_set(PointerRNA *ptr, const char *value)
 	rna_cd_layer_name_set(cd, (CustomDataLayer *)ptr->data, value);
 }
 
+static int rna_Mesh_has_custom_normals_get(PointerRNA *ptr)
+{
+	Mesh *me = ptr->data;
+	return (int)BKE_mesh_has_custom_loop_normals(me);
+}
 
 /* -------------------------------------------------------------------- */
 /* Update Callbacks */
@@ -338,6 +343,17 @@ static void rna_MeshLoop_normal_get(PointerRNA *ptr, float *values)
 	}
 }
 
+static void rna_MeshLoop_normal_set(PointerRNA *ptr, const float *values)
+{
+	Mesh *me = rna_mesh(ptr);
+	MLoop *ml = (MLoop *)ptr->data;
+	float (*vec)[3] = CustomData_get(&me->ldata, (int)(ml - me->mloop), CD_NORMAL);
+
+	if (vec) {
+		normalize_v3_v3(*vec, values);
+	}
+}
+
 static void rna_MeshLoop_tangent_get(PointerRNA *ptr, float *values)
 {
 	Mesh *me = rna_mesh(ptr);
@@ -1974,8 +1990,7 @@ static void rna_def_mloop(BlenderRNA *brna)
 	prop = RNA_def_property(srna, "normal", PROP_FLOAT, PROP_DIRECTION);
 	RNA_def_property_array(prop, 3);
 	RNA_def_property_range(prop, -1.0f, 1.0f);
-	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
-	RNA_def_property_float_funcs(prop, "rna_MeshLoop_normal_get", NULL, NULL);
+	RNA_def_property_float_funcs(prop, "rna_MeshLoop_normal_get", "rna_MeshLoop_normal_set", NULL);
 	RNA_def_property_ui_text(prop, "Normal",
 	                         "Local space unit length split normal vector of this vertex for this polygon "
 	                         "(must be computed beforehand using calc_normals_split or calc_tangents)");
@@ -3186,8 +3201,8 @@ static void rna_def_mesh(BlenderRNA *brna)
 	prop = RNA_def_property(srna, "use_auto_smooth", PROP_BOOLEAN, PROP_NONE);
 	RNA_def_property_boolean_sdna(prop, NULL, "flag", ME_AUTOSMOOTH);
 	RNA_def_property_ui_text(prop, "Auto Smooth",
-	                         "Treat all set-smoothed faces with angles less than the specified angle "
-	                         "as 'smooth', unless they are linked by a sharp edge");
+	                         "Auto smooth (based on smooth/sharp faces/edges and angle between faces), "
+	                         "or use custom split normals data if available");
 	RNA_def_property_update(prop, 0, "rna_Mesh_update_data");
 
 	prop = RNA_def_property(srna, "auto_smooth_angle", PROP_FLOAT, PROP_ANGLE);
@@ -3196,9 +3211,19 @@ static void rna_def_mesh(BlenderRNA *brna)
 	RNA_def_property_range(prop, 0.0f, DEG2RADF(180.0f));
 	RNA_def_property_ui_range(prop, DEG2RADF(0.0f), DEG2RADF(180.0f), 1.0, 1);
 	RNA_def_property_ui_text(prop, "Auto Smooth Angle",
-	                         "Maximum angle between face normals that 'Auto Smooth' will operate on");
+	                         "Maximum angle between face normals that will be considered as smooth "
+	                         "(unused if custom split normals data are available)");
 	RNA_def_property_update(prop, 0, "rna_Mesh_update_data");
 
+	RNA_define_verify_sdna(false);
+	prop = RNA_def_property(srna, "has_custom_normals", PROP_BOOLEAN, PROP_NONE);
+	RNA_def_property_boolean_sdna(prop, NULL, "", 0);
+	RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+	RNA_def_property_ui_text(prop, "Has Custom Normals",
+	                         "True if there are custom split normals data in this mesh");
+	RNA_def_property_boolean_funcs(prop, "rna_Mesh_has_custom_normals_get", NULL);
+	RNA_define_verify_sdna(true);
+
 	prop = RNA_def_property(srna, "show_double_sided", PROP_BOOLEAN, PROP_NONE);
 	RNA_def_property_boolean_sdna(prop, NULL, "flag", ME_TWOSIDED);
 	RNA_def_property_ui_text(prop, "Double Sided", "Render/display the mesh with double or single sided lighting");
diff --git a/source/blender/makesrna/intern/rna_mesh_api.c b/source/blender/makesrna/intern/rna_mesh_api.c
index 09f12f6a4af..40dad902eae 100644
--- a/source/blender/makesrna/intern/rna_mesh_api.c
+++ b/source/blender/makesrna/intern/rna_mesh_api.c
@@ -61,10 +61,19 @@ static const char *rna_Mesh_unit_test_compare(struct Mesh *mesh, struct Mesh *me
 	return ret;
 }
 
+static void rna_Mesh_create_normals_split(Mesh *mesh)
+{
+	if (!CustomData_has_layer(&mesh->ldata, CD_NORMAL)) {
+		CustomData_add_layer(&mesh->ldata, CD_NORMAL, CD_CALLOC, NULL, mesh->totloop);
+		CustomData_set_layer_flag(&mesh->ldata, CD_NORMAL, CD_FLAG_TEMPORARY);
+	}
+}
+
 static void rna_Mesh_calc_normals_split(Mesh *mesh)
 {
 	float (*r_loopnors)[3];
 	float (*polynors)[3];
+	short (*clnors)[2] = NULL;
 	bool free_polynors = false;
 
 	if (CustomData_has_layer(&mesh->ldata, CD_NORMAL)) {
@@ -76,6 +85,9 @@ static void rna_Mesh_calc_normals_split(Mesh *mesh)
 		CustomData_set_layer_flag(&mesh->ldata, CD_NORMAL, CD_FLAG_TEMPORARY);
 	}
 
+	/* may be NULL */
+	clnors = CustomData_get_layer(&mesh->ldata, CD_CUSTOMLOOPNORMAL);
+
 	if (CustomData_has_layer(&mesh->pdata, CD_NORMAL)) {
 		/* This assume that layer is always up to date, not sure this is the case (esp. in Edit mode?)... */
 		polynors = CustomData_get_layer(&mesh->pdata, CD_NORMAL);
@@ -88,9 +100,10 @@ static void rna_Mesh_calc_normals_split(Mesh *mesh)
 		free_polynors = true;
 	}
 
-	BKE_mesh_normals_loop_split(mesh->mvert, mesh->totvert, mesh->medge, mesh->totedge,
-	                            mesh->mloop, r_loopnors, mesh->totloop, mesh->mpoly, polynors, mesh->totpoly,
-	                            (mesh->flag & ME_AUTOSMOOTH) != 0, mesh->smoothresh);
+	BKE_mesh_normals_loop_split(
+	        mesh->mvert, mesh->totvert, mesh->medge, mesh->totedge,
+	        mesh->mloop, r_loopnors, mesh->totloop, mesh->mpoly, (const float (*)[3])polynors, mesh->totpoly,
+	        (mesh->flag & ME_AUTOSMOOTH) != 0, mesh->smoothresh, NULL, clnors, NULL);
 
 	if (free_polynors) {
 		MEM_freeN(polynors);
@@ -139,6 +152,78 @@ static void rna_Mesh_calc_smooth_groups(Mesh *mesh, int use_bitflags, int *r_pol
 	                    r_group_total, use_bitflags);
 }
 
+static void rna_Mesh_normals_split_custom_do(Mesh *mesh, float (*custom_loopnors)[3], const bool use_vertices)
+{
+	float (*polynors)[3];
+	short (*clnors)[2];
+	const int numloops = mesh->totloop;
+	bool free_polynors = false;
+
+	clnors = CustomData_get_layer(&mesh->ldata, CD_CUSTOMLOOPNORMAL);
+	if (clnors) {
+		memset(clnors, 0, sizeof(*clnors) * numloops);
+	}
+	else {
+		clnors = CustomData_add_layer(&mesh->ldata, CD_CUSTOMLOOPNORMAL, CD_DEFAULT, NULL, numloops);
+	}
+
+	if (CustomData_has_layer(&mesh->pdata, CD_NORMAL)) {
+		polynors = CustomData_get_layer(&mesh->pdata, CD_NORMAL);
+	}
+	else {
+		polynors = MEM_mallocN(sizeof(float[3]) * mesh->totpoly, __func__);
+		BKE_mesh_calc_normals_poly(mesh->mvert, mesh->totvert, mesh->mloop, mesh->mpoly, mesh->totloop, mesh->totpoly,
+		                           polynors, false);
+		free_polynors = true;
+	}
+
+	if (use_vertices) {
+		BKE_mesh_normals_loop_custom_from_vertices_set(
+		        mesh->mvert, custom_loopnors, mesh->totvert, mesh->medge, mesh->totedge, mesh->mloop, mesh->totloop,
+		        mesh->mpoly, (const float (*)[3])polynors, mesh->totpoly, clnors);
+	}
+	else {
+		BKE_mesh_normals_loop_custom_set(
+		        mesh->mvert, mesh->totvert, mesh->medge, mesh->totedge, mesh->mloop, custom_loopnors, mesh->totloop,
+		        mesh->mpoly, (const float (*)[3])polynors, mesh->totpoly, clnors);
+	}
+
+	if (free_polynors) {
+		MEM_freeN(polynors);
+	}
+}
+
+static void rna_Mesh_normals_split_custom_set(Mesh *mesh, ReportList *reports, int normals_len, float *normals)
+{
+	float (*loopnors)[3] = (float (*)[3])normals;
+	const int numloops = mesh->totloop;
+
+	if (normals_len != numloops * 3) {
+		BKE_reportf(reports, RPT_ERROR,
+		            "number of custom normals is not number of loops (%f / %d)",
+		            (float)normals_len / 3.0f, numloops);
+		return;
+	}
+
+	rna_Mesh_normals_split_custom_do(mesh, loopnors, false);
+}
+
+static void rna_Mesh_normals_split_custom_set_from_vertices(
+        Mesh *mesh, ReportList *reports, int normals_len, float *normals)
+{
+	float (*vertnors)[3] = (float (*)[3])normals;
+	const int numverts = mesh->totvert;
+
+	if (normals_len != numverts * 3) {
+		BKE_reportf(reports, RPT_ERROR,
+		            "number of custom normals is not number of vertices (%f / %d)",
+		            (float)normals_len / 3.0f, numverts);
+		return;
+	}
+
+	rna_Mesh_normals_split_custom_do(mesh, vertnors, true);
+}
+
 static void rna_Mesh_transform(Mesh *mesh, float *mat, int shape_keys)
 {
 	BKE_mesh_transform(mesh, (float (*)[4])mat, shape_keys);
@@ -152,6 +237,7 @@ void RNA_api_mesh(StructRNA *srna)
 {
 	FunctionRNA *func;
 	PropertyRNA *parm;
+	const int normals_array_dim[] = {1, 3};
 
 	func = RNA_def_function(srna, "transform", "rna_Mesh_transform");
 	RNA_def_function_ui_description(func, "Transform mesh vertices by a matrix");
@@ -162,6 +248,9 @@ void RNA_api_mesh(StructRNA *srna)
 	func = RNA_def_function(srna, "calc_normals", "BKE_mesh_calc_normals");
 	RNA_def_function_ui_description(func, "Calculate vertex normals");
 
+	func = RNA_def_function(srna, "create_normals_split", "rna_Mesh_create_normals_split");
+	RNA_def_function_ui_description(func, "Empty split vertex normals");
+
 	func = RNA_def_function(srna, "calc_normals_split", "rna_Mesh_calc_normals_split");
 	RNA_def_function_ui_description(func, "Calculate split vertex normals, which preserve sharp edges");
 
@@ -192,6 +281,26 @@ void RNA_api_mesh(StructRNA *srna)
 	parm = RNA_def_int(func, "groups", 0, 0, INT_MAX, "groups", "Total number of groups", 0, INT_MAX);
 	RNA_def_property_flag(parm, PROP_OUTPUT);
 
+	func = RNA_def_function(srna, "normals_split_custom_set", "rna_Mesh_normals_split_custom_set");
+	RNA_def_function_ui_description(func,
+	                                "Define custom split normals of this mesh "
+	                                "(use zero-vectors to keep auto ones)");
+	RNA_def_function_flag(func, FUNC_USE_REPORTS);
+	/* TODO, see how array size of 0 works, this shouldnt be used */
+	parm = RNA_def_float_array(func, "normals", 1, NULL, -1.0f, 1.0f, "", "Normals", 0.0f, 0.0f);
+	RNA_def_property_multi_array(parm, 2, normals_array_dim);
+	RNA_def_property_flag(parm, PROP_DYNAMIC | PROP_REQUIRED);
+
+	func = RNA_def_function(srna, "normals_split_custom_set_from_vertices",
+	                        "rna_Mesh_normals_split_custom_set_from_vertices");
+	RNA_def_function_ui_description(func,
+	                                "Define custom split normals of this mesh, from vertices' normals "
+	                                "(use zero-vectors to keep auto ones)");
+	RNA_def_function_flag(func, FUNC_USE_REPORTS);
+	/* TODO, see how array size of 0 works, this shouldnt be used */
+	parm = RNA_def_float_array(func, "normals", 1, NULL, -1.0f, 1.0f, "", "Normals", 0.0f, 0.0f);
+	RNA_def_property_multi_array(parm, 2, normals_array_dim);
+	RNA_def_property_flag(parm, PROP_DYNAMIC | PROP_REQUIRED);
 
 	func = RNA_def_function(srna, "update", "ED_mesh_update");
 	RNA_def_boolean(func, "calc_edges", 0, "Calculate Edges", "Force recalculation of edges");
diff --git a/source/blender/render/intern/source/convertblender.c b/source/blender/render/intern/source/convertblender.c
index 1a8c4681862..ba778ec30e9 100644
--- a/source/blender/render/intern/source/convertblender.c
+++ b/source/blender/render/intern/source/convertblender.c
@@ -3316,7 +3316,7 @@ static void init_render_mesh(Render *re, ObjectRen *obr, int timeoffset)
 							v2= mface->v2;
 							v3= reverse_verts==0 ? mface->v3 : mface->v1;
 							v4= mface->v4;
-							flag= mface->flag & ME_SMOOTH;
+							flag = do_autosmooth ? ME_SMOOTH : mface->flag & ME_SMOOTH;
 
 							vlr= RE_findOrAddVlak(obr, obr->totvlak++);
 							vlr->v1= RE_findOrAddVert(obr, vertofs+v1);