Subsurf: Begin new subdivision surface module

The idea is to use this as a replacement of old CCG, now it is
based on OpenSubdiv. The goal is to reduce any possible overhead
which was happening with OpenSubdiv used by CCG.

Currently implemented/supported:

- Creation from mesh, including topology on OpenSubdiv side,
  its refinement.

- Evaluation of limit point, first order derivatives, normal,
  and face-varying data for individual coarse position.

- Evaluation of whole patches.

  Currently not optimized, uses evaluation of individual coarse
  positions.

- Creation of Mesh from subdiv, with all geometry being real:
  all mvert, medge, mloop, and mpoly.

  This includes custom data interpolation, but all faces currently
  are getting separated (they are converted to ptex patches, which
  we need to weld back together).

Still need to support lighter weights grids and such, but this
is already a required part to have subsurf working in the middle
of modifier stack.

Annoying part is ifdef all over the place, to keep it compilable
when OpenSubdiv is disabled. More cleaner approach would be to
have stub API for OpenSubdiv, so everything gets ifdef-ed in a
much fewer places.

7 files changed, 2120 insertions, 0 deletions

diff --git a/source/blender/blenkernel/BKE_subdiv.h b/source/blender/blenkernel/BKE_subdiv.h
new file mode 100644
index 00000000000..92fb1167f55
--- /dev/null
+++ b/source/blender/blenkernel/BKE_subdiv.h
@@ -0,0 +1,184 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+#ifndef __BKE_SUBDIV_H__
+#define __BKE_SUBDIV_H__
+
+#include "BLI_sys_types.h"
+
+struct Mesh;
+struct OpenSubdiv_Converter;
+struct OpenSubdiv_Evaluator;
+struct OpenSubdiv_TopologyRefiner;
+
+/** \file BKE_subdiv.h
+ *  \ingroup bke
+ *  \since July 2018
+ *  \author Sergey Sharybin
+ */
+
+typedef enum {
+	SUBDIV_FVAR_LINEAR_INTERPOLATION_NONE,
+	SUBDIV_FVAR_LINEAR_INTERPOLATION_CORNERS_ONLY,
+	SUBDIV_FVAR_LINEAR_INTERPOLATION_BOUNDARIES,
+	SUBDIV_FVAR_LINEAR_INTERPOLATION_ALL,
+} eSubdivFVarLinearInterpolation;
+
+typedef struct SubdivSettings {
+	bool is_simple;
+	bool is_adaptive;
+	int level;
+	eSubdivFVarLinearInterpolation fvar_linear_interpolation;
+} SubdivSettings;
+
+typedef struct Subdiv {
+	/* Settings this subdivision surface is created for.
+	 *
+	 * It is read-only after assignment in BKE_subdiv_new_from_FOO().
+	 */
+	SubdivSettings settings;
+
+	/* Total number of ptex faces on subdivision level 0.
+	 *
+	 * Ptex face is what is internally used by OpenSubdiv for evaluator. It is
+	 * a quad face, which corresponds to Blender's legacy Catmull Clark grids.
+	 *
+	 * Basically, here is a correspondence between polygons and ptex faces:
+	 * - Triangle consists of 3 PTex faces.
+	 * - Quad is a single PTex face.
+	 * - N-gon is N PTex faces.
+	 *
+	 * This value is initialized in BKE_subdiv_new_from_FOO() and is read-only
+	 * after this.
+	 */
+	int num_ptex_faces;
+
+	/* Indexed by base face index, element indicates total number of ptex faces
+	 * created for preceding base faces.
+	 */
+	int *face_ptex_offset;
+
+	/* Topology refiner includes all the glue logic to feed Blender side
+	 * topology to OpenSubdiv. It can be shared by both evaluator and GL mesh
+	 * drawer.
+	 */
+	struct OpenSubdiv_TopologyRefiner *topology_refiner;
+
+	/* CPU side evaluator. */
+	struct OpenSubdiv_Evaluator *evaluator;
+} Subdiv;
+
+/* ============================== CONSTRUCTION ============================== */
+
+Subdiv *BKE_subdiv_new_from_converter(const SubdivSettings *settings,
+                                      struct OpenSubdiv_Converter *converter);
+
+Subdiv *BKE_subdiv_new_from_mesh(const SubdivSettings *settings,
+                                 struct Mesh *mesh);
+
+void BKE_subdiv_free(Subdiv *subdiv);
+
+/* ============================= EVALUATION API ============================= */
+
+void BKE_subdiv_eval_begin(Subdiv *subdiv);
+void BKE_subdiv_eval_update_from_mesh(Subdiv *subdiv, const struct Mesh *mesh);
+
+/* Single point queries. */
+
+void BKE_subdiv_eval_limit_point(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3]);
+void BKE_subdiv_eval_limit_point_and_derivatives(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3], float dPdu[3], float dPdv[3]);
+void BKE_subdiv_eval_limit_point_and_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3], float N[3]);
+void BKE_subdiv_eval_limit_point_and_short_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3], short N[3]);
+
+void BKE_subdiv_eval_face_varying(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float varying[2]);
+
+/* Patch queries at given resolution.
+ *
+ * Will evaluate patch at uniformly distributed (u, v) coordinates on a grid
+ * of given resolution, producing resolution^2 evaluation points. The order
+ * goes as u in rows, v in columns.
+ */
+
+void BKE_subdiv_eval_limit_patch_resolution_point(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *buffer, const int offset, const int stride);
+void BKE_subdiv_eval_limit_patch_resolution_point_and_derivatives(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *point_buffer, const int point_offset, const int point_stride,
+        void *du_buffer, const int du_offset, const int du_stride,
+        void *dv_buffer, const int dv_offset, const int dv_stride);
+void BKE_subdiv_eval_limit_patch_resolution_point_and_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *point_buffer, const int point_offset, const int point_stride,
+        void *normal_buffer, const int normal_offset, const int normal_stride);
+void BKE_subdiv_eval_limit_patch_resolution_point_and_short_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *point_buffer, const int point_offset, const int point_stride,
+        void *normal_buffer, const int normal_offset, const int normal_stride);
+
+/* =========================== SUBDIV TO MESH API =========================== */
+
+typedef struct SubdivToMeshSettings {
+	/* Resolution at which ptex are being evaluated.
+	 * This defines how many vertices final mesh will have: every ptex has
+	 * resolution^2 vertices.
+	 */
+	int resolution;
+} SubdivToMeshSettings;
+
+/* Create real hi-res mesh from subdivision, all geometry is "real". */
+struct Mesh *BKE_subdiv_to_mesh(
+        Subdiv *subdiv,
+        const SubdivToMeshSettings *settings,
+        const struct Mesh *coarse_mesh);
+
+#endif  /* __BKE_SUBDIV_H__ */
diff --git a/source/blender/blenkernel/intern/subdiv.c b/source/blender/blenkernel/intern/subdiv.c
new file mode 100644
index 00000000000..72cd39983b9
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv.c
@@ -0,0 +1,114 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/subdiv.c
+ *  \ingroup bke
+ */
+
+#include "BKE_subdiv.h"
+
+#include "BLI_utildefines.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "subdiv_converter.h"
+
+#ifdef WITH_OPENSUBDIV
+#  include "opensubdiv_capi.h"
+#  include "opensubdiv_converter_capi.h"
+#  include "opensubdiv_evaluator_capi.h"
+#  include "opensubdiv_topology_refiner_capi.h"
+#endif
+
+#ifdef WITH_OPENSUBDIV
+static void update_subdiv_after_topology_change(Subdiv *subdiv)
+{
+	/* Count ptex faces. */
+	subdiv->num_ptex_faces = subdiv->topology_refiner->getNumPtexFaces(
+	        subdiv->topology_refiner);
+	/* Initialize offset of base faces in ptex indices. */
+	MEM_SAFE_FREE(subdiv->face_ptex_offset);
+	subdiv->face_ptex_offset = MEM_malloc_arrayN(subdiv->num_ptex_faces,
+	                                             sizeof(int),
+	                                             "subdiv ptex offset");
+	subdiv->topology_refiner->fillFacePtexIndexOffset(
+	        subdiv->topology_refiner,
+	        subdiv->face_ptex_offset);
+}
+#endif
+
+Subdiv *BKE_subdiv_new_from_converter(const SubdivSettings *settings,
+                                      struct OpenSubdiv_Converter *converter)
+{
+#ifdef WITH_OPENSUBDIV
+	OpenSubdiv_TopologyRefinerSettings topology_refiner_settings;
+	topology_refiner_settings.level = settings->level;
+	topology_refiner_settings.is_adaptive = settings->is_adaptive;
+	struct OpenSubdiv_TopologyRefiner *osd_topology_refiner =
+	        openSubdiv_createTopologyRefinerFromConverter(
+	                converter, &topology_refiner_settings);
+	if (osd_topology_refiner == NULL) {
+		return NULL;
+	}
+	Subdiv *subdiv = MEM_callocN(sizeof(Subdiv), "subdiv from converetr");
+	subdiv->settings = *settings;
+	subdiv->topology_refiner = osd_topology_refiner;
+	subdiv->evaluator = NULL;
+	update_subdiv_after_topology_change(subdiv);
+	return subdiv;
+#else
+	UNUSED_VARS(settings, converter);
+	return NULL;
+#endif
+}
+
+Subdiv *BKE_subdiv_new_from_mesh(const SubdivSettings *settings,
+                                 struct Mesh *mesh)
+{
+#ifdef WITH_OPENSUBDIV
+	OpenSubdiv_Converter converter;
+	BKE_subdiv_converter_init_for_mesh(&converter, settings, mesh);
+	Subdiv *subdiv = BKE_subdiv_new_from_converter(settings, &converter);
+	BKE_subdiv_converter_free(&converter);
+	return subdiv;
+#else
+	UNUSED_VARS(settings, mesh);
+	return NULL;
+#endif
+}
+
+void BKE_subdiv_free(Subdiv *subdiv)
+{
+#ifdef WITH_OPENSUBDIV
+	if (subdiv->evaluator != NULL) {
+		openSubdiv_deleteEvaluator(subdiv->evaluator);
+	}
+	if (subdiv->topology_refiner != NULL) {
+		openSubdiv_deleteTopologyRefiner(subdiv->topology_refiner);
+	}
+	MEM_SAFE_FREE(subdiv->face_ptex_offset);
+	MEM_freeN(subdiv);
+#endif
+}
diff --git a/source/blender/blenkernel/intern/subdiv_converter.c b/source/blender/blenkernel/intern/subdiv_converter.c
new file mode 100644
index 00000000000..f6dabfa1c80
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv_converter.c
@@ -0,0 +1,65 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "subdiv_converter.h"
+
+#include "BLI_utildefines.h"
+
+#ifdef WITH_OPENSUBDIV
+#  include "opensubdiv_converter_capi.h"
+#endif
+
+void BKE_subdiv_converter_free(struct OpenSubdiv_Converter *converter)
+{
+#ifdef WITH_OPENSUBDIV
+	if (converter->freeUserData) {
+		converter->freeUserData(converter);
+	}
+#else
+	UNUSED_VARS(converter);
+#endif
+}
+
+/*OpenSubdiv_FVarLinearInterpolation*/ int
+BKE_subdiv_converter_fvar_linear_from_settings(const SubdivSettings *settings)
+{
+#ifdef WITH_OPENSUBDIV
+	switch (settings->fvar_linear_interpolation) {
+		case SUBDIV_FVAR_LINEAR_INTERPOLATION_NONE:
+			return OSD_FVAR_LINEAR_INTERPOLATION_NONE;
+		case SUBDIV_FVAR_LINEAR_INTERPOLATION_CORNERS_ONLY:
+			return OSD_FVAR_LINEAR_INTERPOLATION_CORNERS_ONLY;
+		case SUBDIV_FVAR_LINEAR_INTERPOLATION_BOUNDARIES:
+			return OSD_FVAR_LINEAR_INTERPOLATION_BOUNDARIES;
+		case SUBDIV_FVAR_LINEAR_INTERPOLATION_ALL:
+			return OSD_FVAR_LINEAR_INTERPOLATION_ALL;
+	}
+	BLI_assert(!"Unknown fvar linear interpolation");
+	return OSD_FVAR_LINEAR_INTERPOLATION_NONE;
+#else
+	UNUSED_VARS(settings);
+	return 0;
+#endif
+}
diff --git a/source/blender/blenkernel/intern/subdiv_converter.h b/source/blender/blenkernel/intern/subdiv_converter.h
new file mode 100644
index 00000000000..748d97f4178
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv_converter.h
@@ -0,0 +1,57 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef __BKE_SUBDIV_CONVERTER_H__
+#define __BKE_SUBDIV_CONVERTER_H__
+
+#include "BKE_subdiv.h"
+
+/* NOTE: Was initially used to get proper enumerator types, but this makes
+ * it tricky to compile without OpenSubdiv.
+ */
+/* #include "opensubdiv_converter_capi.h" */
+
+struct Mesh;
+struct OpenSubdiv_Converter;
+struct SubdivSettings;
+
+void BKE_subdiv_converter_init_for_mesh(struct OpenSubdiv_Converter *converter,
+                                        const struct SubdivSettings *settings,
+                                        const struct Mesh *mesh);
+
+/* NOTE: Frees converter data, but not converter itself. This means, that if
+ * converter was allocated on heap, it is up to the user to free that memory.
+ */
+void BKE_subdiv_converter_free(struct OpenSubdiv_Converter *converter);
+
+/* ============================ INTERNAL HELPERS ============================ */
+
+/* TODO(sergey): Find a way to make it OpenSubdiv_FVarLinearInterpolation,
+ * without breaking compilation without OpenSubdiv.
+ */
+int BKE_subdiv_converter_fvar_linear_from_settings(
+        const SubdivSettings *settings);
+
+#endif  /* __BKE_SUBDIV_CONVERTER_H__ */
diff --git a/source/blender/blenkernel/intern/subdiv_converter_mesh.c b/source/blender/blenkernel/intern/subdiv_converter_mesh.c
new file mode 100644
index 00000000000..1a2c26b3564
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv_converter_mesh.c
@@ -0,0 +1,471 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "subdiv_converter.h"
+
+#include <string.h>
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BLI_utildefines.h"
+#include "BLI_math_vector.h"
+
+#include "BKE_customdata.h"
+#include "BKE_mesh_mapping.h"
+#include "BKE_subdiv.h"
+
+#include "MEM_guardedalloc.h"
+
+#ifdef WITH_OPENSUBDIV
+#  include "opensubdiv_capi.h"
+#  include "opensubdiv_converter_capi.h"
+#endif
+
+/* Use mesh element mapping structures during conversion.
+ * Uses more memory but is much faster than naive algorithm.
+ */
+#define USE_MESH_ELEMENT_MAPPING
+
+#ifdef WITH_OPENSUBDIV
+typedef struct ConverterStorage {
+	SubdivSettings settings;
+	const Mesh *mesh;
+
+#ifdef USE_MESH_ELEMENT_MAPPING
+	MeshElemMap *vert_edge_map;
+	MeshElemMap *vert_poly_map;
+	MeshElemMap *edge_poly_map;
+	int *vert_edge_mem;
+	int *vert_poly_mem;
+	int *edge_poly_mem;
+#endif
+
+	/* Indexed by loop index, value denotes index of face-varying vertex
+	 * which corresponds to the UV coordinate.
+	 */
+	int *loop_uv_indices;
+	int num_uv_coordinates;
+} ConverterStorage;
+
+static OpenSubdiv_SchemeType get_scheme_type(
+        const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	if (storage->settings.is_simple) {
+		return OSD_SCHEME_BILINEAR;
+	}
+	else {
+		return OSD_SCHEME_CATMARK;
+	}
+}
+
+static OpenSubdiv_FVarLinearInterpolation get_fvar_linear_interpolation(
+        const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	return BKE_subdiv_converter_fvar_linear_from_settings(&storage->settings);
+}
+
+static int get_num_faces(const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	return storage->mesh->totpoly;
+}
+
+static int get_num_edges(const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	return storage->mesh->totedge;
+}
+
+static int get_num_verts(const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	return storage->mesh->totvert;
+}
+
+static int get_num_face_verts(const OpenSubdiv_Converter *converter, int face)
+{
+	ConverterStorage *storage = converter->user_data;
+	return storage->mesh->mpoly[face].totloop;
+}
+
+static void get_face_verts(const OpenSubdiv_Converter *converter,
+                           int face,
+                           int *face_verts)
+{
+	ConverterStorage *storage = converter->user_data;
+	const MPoly *mp = &storage->mesh->mpoly[face];
+	const MLoop *mloop = storage->mesh->mloop;
+	for (int loop = 0; loop < mp->totloop; loop++) {
+		face_verts[loop] = mloop[mp->loopstart + loop].v;
+	}
+}
+
+static void get_face_edges(const OpenSubdiv_Converter *converter,
+                           int face,
+                           int *face_edges)
+{
+	ConverterStorage *storage = converter->user_data;
+	const MPoly *mp = &storage->mesh->mpoly[face];
+	const MLoop *mloop = storage->mesh->mloop;
+	for (int loop = 0; loop < mp->totloop; loop++) {
+		face_edges[loop] = mloop[mp->loopstart + loop].e;
+	}
+}
+
+static void get_edge_verts(const OpenSubdiv_Converter *converter,
+                           int edge,
+                           int *edge_verts)
+{
+	ConverterStorage *storage = converter->user_data;
+	const MEdge *me = &storage->mesh->medge[edge];
+	edge_verts[0] = me->v1;
+	edge_verts[1] = me->v2;
+}
+
+static int get_num_edge_faces(const OpenSubdiv_Converter *converter, int edge)
+{
+	ConverterStorage *storage = converter->user_data;
+#ifdef USE_MESH_ELEMENT_MAPPING
+	return storage->edge_poly_map[edge].count;
+#else
+	const Mesh *mesh = storage->mesh;
+	const MPoly *mpoly = mesh->mpoly;
+	const MLoop *mloop = mesh->mloop;
+	int num = 0;
+	for (int poly = 0; poly < mesh->totpoly; poly++) {
+		const MPoly *mp = &mpoly[poly];
+		for (int loop = 0; loop < mp->totloop; loop++) {
+			const MLoop *ml = &mloop[mp->loopstart + loop];
+			if (ml->e == edge) {
+				++num;
+				break;
+			}
+		}
+	}
+	return num;
+#endif
+}
+
+static void get_edge_faces(const OpenSubdiv_Converter *converter,
+                           int edge,
+                           int *edge_faces)
+{
+	ConverterStorage *storage = converter->user_data;
+#ifdef USE_MESH_ELEMENT_MAPPING
+	memcpy(edge_faces,
+	       storage->edge_poly_map[edge].indices,
+	       sizeof(int) * storage->edge_poly_map[edge].count);
+#else
+	const Mesh *mesh = storage->mesh;
+	const MPoly *mpoly = mesh->mpoly;
+	const MLoop *mloop = mesh->mloop;
+	int num = 0;
+	for (int poly = 0; poly < mesh->totpoly; poly++) {
+		const MPoly *mp = &mpoly[poly];
+		for (int loop = 0; loop < mpoly->totloop; loop++) {
+			const MLoop *ml = &mloop[mp->loopstart + loop];
+			if (ml->e == edge) {
+				edge_faces[num++] = poly;
+				break;
+			}
+		}
+	}
+#endif
+}
+
+static float get_edge_sharpness(const OpenSubdiv_Converter *converter, int edge)
+{
+	ConverterStorage *storage = converter->user_data;
+	const MEdge *medge = storage->mesh->medge;
+	const float edge_crease =  (float)medge[edge].crease / 255.0f;
+	return edge_crease * storage->settings.level;
+}
+
+static int get_num_vert_edges(const OpenSubdiv_Converter *converter, int vert)
+{
+	ConverterStorage *storage = converter->user_data;
+#ifdef USE_MESH_ELEMENT_MAPPING
+	return storage->vert_edge_map[vert].count;
+#else
+	const Mesh *mesh = storage->mesh;
+	const MEdge *medge = mesh->medge;
+	int num = 0;
+	for (int edge = 0; edge < mesh->totedge; edge++) {
+		const MEdge *me = &medge[edge];
+		if (me->v1 == vert || me->v2 == vert) {
+			++num;
+		}
+	}
+	return num;
+#endif
+}
+
+static void get_vert_edges(const OpenSubdiv_Converter *converter,
+                           int vert,
+                           int *vert_edges)
+{
+	ConverterStorage *storage = converter->user_data;
+#ifdef USE_MESH_ELEMENT_MAPPING
+	memcpy(vert_edges,
+	       storage->vert_edge_map[vert].indices,
+	       sizeof(int) * storage->vert_edge_map[vert].count);
+#else
+	const Mesh *mesh = storage->mesh;
+	const MEdge *medge = mesh->medge;
+	int num = 0;
+	for (int edge = 0; edge < mesh->totedge; edge++) {
+		const MEdge *me = &medge[edge];
+		if (me->v1 == vert || me->v2 == vert) {
+			vert_edges[num++] = edge;
+		}
+	}
+#endif
+}
+
+static int get_num_vert_faces(const OpenSubdiv_Converter *converter, int vert)
+{
+	ConverterStorage *storage = converter->user_data;
+#ifdef USE_MESH_ELEMENT_MAPPING
+	return storage->vert_poly_map[vert].count;
+#else
+	const Mesh *mesh = storage->mesh;
+	const MPoly *mpoly = mesh->mpoly;
+	const MLoop *mloop = mesh->mloop;
+	int num = 0;
+	for (int poly = 0; poly < mesh->totpoly; poly++) {
+		const MPoly *mp = &mpoly[poly];
+		for (int loop = 0; loop < mpoly->totloop; loop++) {
+			const MLoop *ml = &mloop[mp->loopstart + loop];
+			if (ml->v == vert) {
+				++num;
+				break;
+			}
+		}
+	}
+	return num;
+#endif
+}
+
+static void get_vert_faces(const OpenSubdiv_Converter *converter,
+                           int vert,
+                           int *vert_faces)
+{
+	ConverterStorage *storage = converter->user_data;
+#ifdef USE_MESH_ELEMENT_MAPPING
+	memcpy(vert_faces,
+	       storage->vert_poly_map[vert].indices,
+	       sizeof(int) * storage->vert_poly_map[vert].count);
+#else
+	const Mesh *mesh = storage->mesh;
+	const MPoly *mpoly = mesh->mpoly;
+	const MLoop *mloop = mesh->mloop;
+	int num = 0;
+	for (int poly = 0; poly < mesh->totpoly; poly++) {
+		const MPoly *mp = &mpoly[poly];
+		for (int loop = 0; loop < mpoly->totloop; loop++) {
+			const MLoop *ml = &mloop[mp->loopstart + loop];
+			if (ml->v == vert) {
+				vert_faces[num++] = poly;
+				break;
+			}
+		}
+	}
+#endif
+}
+
+static int get_num_uv_layers(const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	const Mesh *mesh = storage->mesh;
+	return CustomData_number_of_layers(&mesh->ldata, CD_MLOOPUV);
+}
+
+static void precalc_uv_layer(const OpenSubdiv_Converter *converter,
+                             const int layer_index)
+{
+	ConverterStorage *storage = converter->user_data;
+	const Mesh *mesh = storage->mesh;
+	const MPoly *mpoly = mesh->mpoly;
+	const MLoop *mloop = mesh->mloop;
+	const MLoopUV *mloopuv = CustomData_get_layer_n(
+	        &mesh->ldata, CD_MLOOPUV, layer_index);
+	const int num_poly = mesh->totpoly;
+	const int num_vert = mesh->totvert;
+	const float limit[2] = {STD_UV_CONNECT_LIMIT, STD_UV_CONNECT_LIMIT};
+	/* Initialize memory required for the operations. */
+	if (storage->loop_uv_indices == NULL) {
+		storage->loop_uv_indices = MEM_malloc_arrayN(
+		        mesh->totloop, sizeof(int), "loop uv vertex index");
+	}
+	UvVertMap *uv_vert_map = BKE_mesh_uv_vert_map_create(
+	        mpoly, mloop, mloopuv,
+	        num_poly, num_vert,
+	        limit,
+	        false, true);
+	/* NOTE: First UV vertex is supposed to be always marked as separate. */
+	storage->num_uv_coordinates = -1;
+	for (int vertex_index = 0; vertex_index < num_vert; ++vertex_index) {
+		const UvMapVert *uv_vert = BKE_mesh_uv_vert_map_get_vert(uv_vert_map,
+		                                                         vertex_index);
+		while (uv_vert != NULL) {
+			if (uv_vert->separate) {
+				storage->num_uv_coordinates++;
+			}
+			const MPoly *mp = &mpoly[uv_vert->poly_index];
+			const int global_loop_index = mp->loopstart +
+			                              uv_vert->loop_of_poly_index;
+			storage->loop_uv_indices[global_loop_index] =
+			        storage->num_uv_coordinates;
+			uv_vert = uv_vert->next;
+		}
+	}
+	/* So far this value was used as a 0-based index, actual number of UV
+	 * vertices is 1 more.
+	 */
+	storage->num_uv_coordinates += 1;
+	BKE_mesh_uv_vert_map_free(uv_vert_map);
+}
+
+static void finish_uv_layer(const OpenSubdiv_Converter *UNUSED(converter))
+{
+}
+
+static int get_num_uvs(const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *storage = converter->user_data;
+	return storage->num_uv_coordinates;
+}
+
+static int get_face_corner_uv_index(const OpenSubdiv_Converter *converter,
+                                    const int face_index,
+                                    const int corner)
+{
+	ConverterStorage *storage = converter->user_data;
+	const MPoly *mp = &storage->mesh->mpoly[face_index];
+	return storage->loop_uv_indices[mp->loopstart + corner];
+}
+
+static void free_user_data(const OpenSubdiv_Converter *converter)
+{
+	ConverterStorage *user_data = converter->user_data;
+	MEM_SAFE_FREE(user_data->loop_uv_indices);
+#ifdef USE_MESH_ELEMENT_MAPPING
+	MEM_freeN(user_data->vert_edge_map);
+	MEM_freeN(user_data->vert_edge_mem);
+	MEM_freeN(user_data->vert_poly_map);
+	MEM_freeN(user_data->vert_poly_mem);
+	MEM_freeN(user_data->edge_poly_map);
+	MEM_freeN(user_data->edge_poly_mem);
+#endif
+	MEM_freeN(user_data);
+}
+
+static void init_functions(OpenSubdiv_Converter *converter)
+{
+	converter->getSchemeType = get_scheme_type;
+
+	converter->getFVarLinearInterpolation = get_fvar_linear_interpolation;
+
+	converter->getNumFaces = get_num_faces;
+	converter->getNumEdges = get_num_edges;
+	converter->getNumVertices = get_num_verts;
+
+	converter->getNumFaceVertices = get_num_face_verts;
+	converter->getFaceVertices = get_face_verts;
+	converter->getFaceEdges = get_face_edges;
+
+	converter->getEdgeVertices = get_edge_verts;
+	converter->getNumEdgeFaces = get_num_edge_faces;
+	converter->getEdgeFaces = get_edge_faces;
+	converter->getEdgeSharpness = get_edge_sharpness;
+
+	converter->getNumVertexEdges = get_num_vert_edges;
+	converter->getVertexEdges = get_vert_edges;
+	converter->getNumVertexFaces = get_num_vert_faces;
+	converter->getVertexFaces = get_vert_faces;
+
+	converter->getNumUVLayers = get_num_uv_layers;
+	converter->precalcUVLayer = precalc_uv_layer;
+	converter->finishUVLayer = finish_uv_layer;
+	converter->getNumUVCoordinates = get_num_uvs;
+	converter->getFaceCornerUVIndex = get_face_corner_uv_index;
+
+	converter->freeUserData = free_user_data;
+}
+
+static void create_element_maps_if_needed(ConverterStorage *storage)
+{
+#ifdef USE_MESH_ELEMENT_MAPPING
+	const Mesh *mesh = storage->mesh;
+	BKE_mesh_vert_edge_map_create(&storage->vert_edge_map,
+	                              &storage->vert_edge_mem,
+	                              mesh->medge,
+	                              mesh->totvert,
+	                              mesh->totedge);
+	BKE_mesh_vert_poly_map_create(&storage->vert_poly_map,
+	                              &storage->vert_poly_mem,
+	                              mesh->mpoly,
+	                              mesh->mloop,
+	                              mesh->totvert,
+	                              mesh->totpoly,
+	                              mesh->totloop);
+	BKE_mesh_edge_poly_map_create(&storage->edge_poly_map,
+	                              &storage->edge_poly_mem,
+	                              mesh->medge, mesh->totedge,
+	                              mesh->mpoly, mesh->totpoly,
+	                              mesh->mloop, mesh->totloop);
+#else
+	(void) storage;  /* Ignored. */
+#endif
+}
+
+static void init_user_data(OpenSubdiv_Converter *converter,
+                           const SubdivSettings *settings,
+                           const Mesh *mesh)
+{
+	ConverterStorage *user_data =
+	         MEM_mallocN(sizeof(ConverterStorage), __func__);
+	user_data->settings = *settings;
+	user_data->mesh = mesh;
+	user_data->loop_uv_indices = NULL;
+	create_element_maps_if_needed(user_data);
+	converter->user_data = user_data;
+}
+#endif
+
+void BKE_subdiv_converter_init_for_mesh(struct OpenSubdiv_Converter *converter,
+                                        const SubdivSettings *settings,
+                                        const Mesh *mesh)
+{
+#ifdef WITH_OPENSUBDIV
+	init_functions(converter);
+	init_user_data(converter, settings, mesh);
+#else
+	UNUSED_VARS(converter, settings, mesh);
+#endif
+}
diff --git a/source/blender/blenkernel/intern/subdiv_eval.c b/source/blender/blenkernel/intern/subdiv_eval.c
new file mode 100644
index 00000000000..0ab19fc8df5
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv_eval.c
@@ -0,0 +1,319 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/subdiv_eval.c
+ *  \ingroup bke
+ */
+
+#include "BKE_subdiv.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BLI_utildefines.h"
+#include "BLI_math_vector.h"
+
+#include "BKE_customdata.h"
+
+#ifdef WITH_OPENSUBDIV
+#  include "opensubdiv_evaluator_capi.h"
+#  include "opensubdiv_topology_refiner_capi.h"
+#endif
+
+void BKE_subdiv_eval_begin(Subdiv *subdiv)
+{
+#ifdef WITH_OPENSUBDIV
+	if (subdiv->evaluator == NULL) {
+		subdiv->evaluator = openSubdiv_createEvaluatorFromTopologyRefiner(
+		        subdiv->topology_refiner);
+	}
+	else {
+		/* TODO(sergey): Check for topology change. */
+	}
+#else
+	UNUSED_VARS(subdiv);
+#endif
+}
+
+#ifdef WITH_OPENSUBDIV
+static void set_face_varying_data_from_uv(Subdiv *subdiv,
+                                          const MLoopUV *mloopuv,
+                                          const int layer_index)
+{
+	OpenSubdiv_TopologyRefiner *topology_refiner = subdiv->topology_refiner;
+	OpenSubdiv_Evaluator *evaluator = subdiv->evaluator;
+	const int num_faces = topology_refiner->getNumFaces(topology_refiner);
+	const MLoopUV *mluv = mloopuv;
+	/* TODO(sergey): OpenSubdiv's C-API converter can change winding of
+	 * loops of a face, need to watch for that, to prevent wrong UVs assigned.
+	 */
+	for (int face_index = 0; face_index < num_faces; ++face_index) {
+		const int num_face_vertices = topology_refiner->getNumFaceVertices(
+		        topology_refiner, face_index);
+		const int *uv_indicies = topology_refiner->getFaceFVarValueIndices(
+		        topology_refiner, face_index, layer_index);
+		for (int vertex_index = 0;
+		     vertex_index < num_face_vertices;
+		     vertex_index++, mluv++)
+		{
+			evaluator->setFaceVaryingData(evaluator,
+			                              mluv->uv,
+			                              uv_indicies[vertex_index],
+			                              1);
+		}
+	}
+}
+#endif
+
+void BKE_subdiv_eval_update_from_mesh(Subdiv *subdiv, const Mesh *mesh)
+{
+#ifdef WITH_OPENSUBDIV
+	BKE_subdiv_eval_begin(subdiv);
+	/* Set coordinates of base mesh vertices. */
+	subdiv->evaluator->setCoarsePositionsFromBuffer(
+	        subdiv->evaluator,
+	        mesh->mvert,
+	        offsetof(MVert, co),
+	        sizeof(MVert),
+	        0, mesh->totvert);
+	/* Set face-varyign data to UV maps. */
+	const int num_uv_layers =
+	        CustomData_number_of_layers(&mesh->ldata, CD_MLOOPUV);
+	for (int layer_index = 0; layer_index < num_uv_layers; layer_index++) {
+		const MLoopUV *mloopuv = CustomData_get_layer_n(
+		        &mesh->ldata, CD_MLOOPUV, layer_index);
+		set_face_varying_data_from_uv(subdiv, mloopuv, layer_index);
+		/* NOTE: Currently evaluator can only handle single face varying layer.
+		 * This is a limitation of C-API and some underlying helper classes from
+		 * our side which will get fixed.
+		 */
+		break;
+	}
+	/* Update evaluator to the new coarse geometry. */
+	subdiv->evaluator->refine(subdiv->evaluator);
+#else
+	UNUSED_VARS(subdiv, mesh);
+#endif
+}
+
+/* ========================== Single point queries ========================== */
+
+void BKE_subdiv_eval_limit_point(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3])
+{
+	BKE_subdiv_eval_limit_point_and_derivatives(subdiv,
+	                                            ptex_face_index,
+	                                            u, v,
+	                                            P, NULL, NULL);
+}
+
+void BKE_subdiv_eval_limit_point_and_derivatives(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3], float dPdu[3], float dPdv[3])
+{
+#ifdef WITH_OPENSUBDIV
+	subdiv->evaluator->evaluateLimit(subdiv->evaluator,
+	                                 ptex_face_index,
+	                                 u, v,
+	                                 P, dPdu, dPdv);
+#else
+	UNUSED_VARS(subdiv, ptex_face_index, u, v, P, dPdu, dPdv);
+#endif
+}
+
+void BKE_subdiv_eval_limit_point_and_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3], float N[3])
+{
+	float dPdu[3], dPdv[3];
+	BKE_subdiv_eval_limit_point_and_derivatives(subdiv,
+	                                            ptex_face_index,
+	                                            u, v,
+	                                            P, dPdu, dPdv);
+	cross_v3_v3v3(N, dPdu, dPdv);
+	normalize_v3(N);
+}
+
+void BKE_subdiv_eval_limit_point_and_short_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float P[3], short N[3])
+{
+	float N_float[3];
+	BKE_subdiv_eval_limit_point_and_normal(subdiv,
+	                                       ptex_face_index,
+	                                       u, v,
+	                                       P, N_float);
+	normal_float_to_short_v3(N, N_float);
+}
+
+void BKE_subdiv_eval_face_varying(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const float u, const float v,
+        float face_varying[2])
+{
+#ifdef WITH_OPENSUBDIV
+	subdiv->evaluator->evaluateFaceVarying(subdiv->evaluator,
+	                                       ptex_face_index,
+	                                       u, v,
+	                                       face_varying);
+#else
+	UNUSED_VARS(subdiv, ptex_face_index, u, v, face_varying);
+#endif
+}
+
+/* ===================  Patch queries at given resolution =================== */
+
+/* Move buffer forward by a given number of bytes. */
+static void buffer_apply_offset(void **buffer, const int offset)
+{
+  *buffer = ((unsigned char *)*buffer) + offset;
+}
+
+/* Write given number of floats to the beginning of given buffer.  */
+static void buffer_write_float_value(void **buffer,
+                                     const float *values_buffer, int num_values)
+{
+	memcpy(*buffer, values_buffer, sizeof(float) * num_values);
+}
+
+/* Similar to above, just operates with short values. */
+static void buffer_write_short_value(void **buffer,
+                                     const short *values_buffer, int num_values)
+{
+	memcpy(*buffer, values_buffer, sizeof(short) * num_values);
+}
+
+void BKE_subdiv_eval_limit_patch_resolution_point(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *buffer, const int offset, const int stride)
+{
+	buffer_apply_offset(&buffer, offset);
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	for (int y = 0; y < resolution; y++) {
+		const float v = y * inv_resolution_1;
+		for (int x = 0; x < resolution; x++) {
+			const float u = x * inv_resolution_1;
+			BKE_subdiv_eval_limit_point(subdiv,
+			                            ptex_face_index,
+			                            u, v,
+			                            buffer);
+			buffer_apply_offset(&buffer, stride);
+		}
+	}
+}
+
+void BKE_subdiv_eval_limit_patch_resolution_point_and_derivatives(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *point_buffer, const int point_offset, const int point_stride,
+        void *du_buffer, const int du_offset, const int du_stride,
+        void *dv_buffer, const int dv_offset, const int dv_stride)
+{
+	buffer_apply_offset(&point_buffer, point_offset);
+	buffer_apply_offset(&du_buffer, du_offset);
+	buffer_apply_offset(&dv_buffer, dv_offset);
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	for (int y = 0; y < resolution; y++) {
+		const float v = y * inv_resolution_1;
+		for (int x = 0; x < resolution; x++) {
+			const float u = x * inv_resolution_1;
+			BKE_subdiv_eval_limit_point_and_derivatives(
+			        subdiv,
+			        ptex_face_index,
+			        u, v,
+			        point_buffer, du_buffer, dv_buffer);
+			buffer_apply_offset(&point_buffer, point_stride);
+			buffer_apply_offset(&du_buffer, du_stride);
+			buffer_apply_offset(&dv_buffer, dv_stride);
+		}
+	}
+}
+
+void BKE_subdiv_eval_limit_patch_resolution_point_and_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *point_buffer, const int point_offset, const int point_stride,
+        void *normal_buffer, const int normal_offset, const int normal_stride)
+{
+	buffer_apply_offset(&point_buffer, point_offset);
+	buffer_apply_offset(&normal_buffer, normal_offset);
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	for (int y = 0; y < resolution; y++) {
+		const float v = y * inv_resolution_1;
+		for (int x = 0; x < resolution; x++) {
+			const float u = x * inv_resolution_1;
+			float normal[3];
+			BKE_subdiv_eval_limit_point_and_normal(
+			        subdiv,
+			        ptex_face_index,
+			        u, v,
+			        point_buffer, normal);
+			buffer_write_float_value(&normal_buffer, normal, 3);
+			buffer_apply_offset(&point_buffer, point_stride);
+			buffer_apply_offset(&normal_buffer, normal_stride);
+		}
+	}
+}
+
+void BKE_subdiv_eval_limit_patch_resolution_point_and_short_normal(
+        Subdiv *subdiv,
+        const int ptex_face_index,
+        const int resolution,
+        void *point_buffer, const int point_offset, const int point_stride,
+        void *normal_buffer, const int normal_offset, const int normal_stride)
+{
+	buffer_apply_offset(&point_buffer, point_offset);
+	buffer_apply_offset(&normal_buffer, normal_offset);
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	for (int y = 0; y < resolution; y++) {
+		const float v = y * inv_resolution_1;
+		for (int x = 0; x < resolution; x++) {
+			const float u = x * inv_resolution_1;
+			short normal[3];
+			BKE_subdiv_eval_limit_point_and_short_normal(
+			        subdiv,
+			        ptex_face_index,
+			        u, v,
+			        point_buffer, normal);
+			buffer_write_short_value(&normal_buffer, normal, 3);
+			buffer_apply_offset(&point_buffer, point_stride);
+			buffer_apply_offset(&normal_buffer, normal_stride);
+		}
+	}
+}
diff --git a/source/blender/blenkernel/intern/subdiv_mesh.c b/source/blender/blenkernel/intern/subdiv_mesh.c
new file mode 100644
index 00000000000..8a58605ba19
--- /dev/null
+++ b/source/blender/blenkernel/intern/subdiv_mesh.c
@@ -0,0 +1,910 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2018 by Blender Foundation.
+ * All rights reserved.
+ *
+ * Contributor(s): Sergey Sharybin.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/subdiv_mesh.c
+ *  \ingroup bke
+ */
+
+#include "BKE_subdiv.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "BLI_alloca.h"
+#include "BLI_math_vector.h"
+#include "BLI_task.h"
+
+#include "BKE_mesh.h"
+
+/* TODO(sergey): Somehow move this to subdiv code? */
+static int mpoly_ptex_faces_count_get(const MPoly *mp)
+{
+	if (mp->totloop == 4) {
+		return 1;
+	}
+	else {
+		return mp->totloop;
+	}
+}
+
+static int num_edges_per_ptex_get(const int resolution)
+{
+	return 2 * (resolution - 1) * resolution;
+}
+
+static int num_polys_per_ptex_get(const int resolution)
+{
+	return (resolution - 1) * (resolution - 1);
+}
+
+typedef struct SubdivMeshContext {
+	const Mesh *coarse_mesh;
+	Subdiv *subdiv;
+	Mesh *subdiv_mesh;
+	const SubdivToMeshSettings *settings;
+	/* Cached custom data arrays for fastter access. */
+	int *vert_origindex;
+	int *edge_origindex;
+	int *loop_origindex;
+	int *poly_origindex;
+	/* UV layers interpolation. */
+	int num_uv_layers;
+	MLoopUV *uv_layers[MAX_MTFACE];
+} SubdivMeshContext;
+
+typedef struct LoopsOfPtex {
+	/* First loop of the ptex, starts at ptex (0, 0) and goes in u direction. */
+	const MLoop *first_loop;
+	/* Last loop of the ptex, starts at ptex (0, 0) and goes in v direction. */
+	const MLoop *last_loop;
+	/* For quad coarse faces only. */
+	const MLoop *second_loop;
+	const MLoop *third_loop;
+} LoopsOfPtex;
+
+static void loops_of_ptex_get(
+        const SubdivMeshContext *ctx,
+        LoopsOfPtex *loops_of_ptex,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	const MLoop *coarse_mloop = ctx->coarse_mesh->mloop;
+	const int first_ptex_loop_index = coarse_poly->loopstart + ptex_face_index;
+	/* Loop which look in the (opposite) V direction of the current
+	 * ptex face.
+	 *
+	 * TOOD(sergey): Get rid of using module on every iteration.
+	 */
+	const int last_ptex_loop_index =
+	        coarse_poly->loopstart +
+	        (ptex_face_index + coarse_poly->totloop - 1) % coarse_poly->totloop;
+	loops_of_ptex->first_loop = &coarse_mloop[first_ptex_loop_index];
+	loops_of_ptex->last_loop = &coarse_mloop[last_ptex_loop_index];
+	if (coarse_poly->totloop == 4) {
+		loops_of_ptex->second_loop = loops_of_ptex->first_loop + 1;
+		loops_of_ptex->third_loop = loops_of_ptex->first_loop + 2;
+	}
+	else {
+		loops_of_ptex->second_loop = NULL;
+		loops_of_ptex->third_loop = NULL;
+	}
+}
+
+typedef struct EdgesOfPtex {
+	/* First edge of the ptex, starts at ptex (0, 0) and goes in u direction. */
+	const MEdge *first_edge;
+	/* Last edge of the ptex, starts at ptex (0, 0) and goes in v direction. */
+	const MEdge *last_edge;
+	/* For quad coarse faces only. */
+	const MEdge *second_edge;
+	const MEdge *third_edge;
+} EdgesOfPtex;
+
+static void edges_of_ptex_get(
+        const SubdivMeshContext *ctx,
+        EdgesOfPtex *edges_of_ptex,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	const MEdge *coarse_medge = ctx->coarse_mesh->medge;
+	LoopsOfPtex loops_of_ptex;
+	loops_of_ptex_get(ctx, &loops_of_ptex, coarse_poly, ptex_face_index);
+	edges_of_ptex->first_edge = &coarse_medge[loops_of_ptex.first_loop->e];
+	edges_of_ptex->last_edge = &coarse_medge[loops_of_ptex.last_loop->e];
+	if (coarse_poly->totloop == 4) {
+		edges_of_ptex->second_edge =
+		        &coarse_medge[loops_of_ptex.second_loop->e];
+		edges_of_ptex->third_edge =
+		        &coarse_medge[loops_of_ptex.third_loop->e];
+	}
+	else {
+		edges_of_ptex->second_edge = NULL;
+		edges_of_ptex->third_edge = NULL;
+	}
+}
+
+/* TODO(sergey): Somehow de-duplicate with loops storage, without too much
+ * exception cases all over the code.
+ */
+
+typedef struct VerticesForInterpolation {
+	/* This field points to a vertex data which is to be used for interpolation.
+	 * The idea is to avoid unnecessary allocations for regular faces, where
+	 * we can simply 
+	 */
+	const CustomData *vertex_data;
+	/* Vertices data calculated for ptex corners. There are always 4 elements
+	 * in this custom data, aligned the following way:
+	 *
+	 *   index 0 -> uv (0, 0)
+	 *   index 1 -> uv (0, 1)
+	 *   index 2 -> uv (1, 1)
+	 *   index 3 -> uv (1, 0)
+	 *
+	 * Is allocated for non-regular faces (triangles and n-gons).
+	 */
+	CustomData vertex_data_storage;
+	bool vertex_data_storage_allocated;
+	/* Infices within vertex_data to interpolate for. The indices are aligned
+	 * with uv coordinates in a similar way as indices in loop_data_storage.
+	 */
+	int vertex_indices[4];
+} VerticesForInterpolation;
+
+static void vertex_interpolation_init(
+        const SubdivMeshContext *ctx,
+        VerticesForInterpolation *vertex_interpolation,
+		const MPoly *coarse_poly)
+{
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MLoop *coarse_mloop = coarse_mesh->mloop;
+	if (coarse_poly->totloop == 4) {
+		vertex_interpolation->vertex_data = &coarse_mesh->vdata;
+		vertex_interpolation->vertex_indices[0] =
+		        coarse_mloop[coarse_poly->loopstart + 0].v;
+		vertex_interpolation->vertex_indices[1] =
+		        coarse_mloop[coarse_poly->loopstart + 1].v;
+		vertex_interpolation->vertex_indices[2] =
+		        coarse_mloop[coarse_poly->loopstart + 2].v;
+		vertex_interpolation->vertex_indices[3] =
+		        coarse_mloop[coarse_poly->loopstart + 3].v;
+		vertex_interpolation->vertex_data_storage_allocated = false;
+	}
+	else {
+		vertex_interpolation->vertex_data =
+		        &vertex_interpolation->vertex_data_storage;
+		/* Allocate storage for loops corresponding to ptex corners. */
+		CustomData_copy(&ctx->coarse_mesh->vdata,
+		                &vertex_interpolation->vertex_data_storage,
+		                CD_MASK_EVERYTHING,
+		                CD_CALLOC,
+		                4);
+		/* Initialize indices. */
+		vertex_interpolation->vertex_indices[0] = 0;
+		vertex_interpolation->vertex_indices[1] = 1;
+		vertex_interpolation->vertex_indices[2] = 2;
+		vertex_interpolation->vertex_indices[3] = 3;
+		vertex_interpolation->vertex_data_storage_allocated = true;
+		/* Interpolate center of poly right away, it stays unchanged for all
+		 * ptex faces.
+		 */
+		const float weight = 1.0f / (float)coarse_poly->totloop;
+		float *weights = BLI_array_alloca(weights, coarse_poly->totloop);
+		int *indices = BLI_array_alloca(indices, coarse_poly->totloop);
+		for (int i = 0; i < coarse_poly->totloop; ++i) {
+			weights[i] = weight;
+			indices[i] = coarse_poly->loopstart + i;
+		}
+		CustomData_interp(&coarse_mesh->vdata,
+		                  &vertex_interpolation->vertex_data_storage,
+		                  indices,
+		                  weights, NULL,
+		                  coarse_poly->totloop,
+		                  2);
+	}
+}
+
+static void vertex_interpolation_from_ptex(
+        const SubdivMeshContext *ctx,
+        VerticesForInterpolation *vertex_interpolation,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	if (coarse_poly->totloop == 4) {
+		/* Nothing to do, all indices and data is already assigned. */
+	} else {
+		const CustomData *vertex_data = &ctx->coarse_mesh->vdata;
+		const Mesh *coarse_mesh = ctx->coarse_mesh;
+		const MLoop *coarse_mloop = coarse_mesh->mloop;
+		LoopsOfPtex loops_of_ptex;
+		loops_of_ptex_get(ctx, &loops_of_ptex, coarse_poly, ptex_face_index);
+		/* Ptex face corner corresponds to a poly loop with same index. */
+		CustomData_copy_data(
+		        vertex_data,
+		        &vertex_interpolation->vertex_data_storage,
+		        coarse_mloop[coarse_poly->loopstart + ptex_face_index].v,
+		        0,
+		        1);
+		/* Interpolate remaining ptex face corners, which hits loops
+		 * middle points.
+		 *
+		 * TODO(sergey): Re-use one of interpolation results from previous
+		 * iteration.
+		 */
+		const float weights[2] = {0.5f, 0.5f};
+		const int first_indices[2] = {
+		        coarse_mloop[loops_of_ptex.first_loop - coarse_mloop].v,
+		        coarse_mloop[(loops_of_ptex.first_loop + 1 - coarse_mloop) %
+		                coarse_poly->totloop].v};
+		const int last_indices[2] = {
+		        coarse_mloop[loops_of_ptex.last_loop - coarse_mloop].v,
+		        coarse_mloop[loops_of_ptex.first_loop - coarse_mloop].v};
+		CustomData_interp(vertex_data,
+		                  &vertex_interpolation->vertex_data_storage,
+		                  first_indices,
+		                  weights, NULL,
+		                  2,
+		                  1);
+		CustomData_interp(vertex_data,
+		                  &vertex_interpolation->vertex_data_storage,
+		                  last_indices,
+		                  weights, NULL,
+		                  2,
+		                  3);
+	}
+}
+
+static void vertex_interpolation_end(
+        VerticesForInterpolation *vertex_interpolation)
+{
+	if (vertex_interpolation->vertex_data_storage_allocated) {
+		CustomData_free(&vertex_interpolation->vertex_data_storage, 4);
+	}
+}
+
+typedef struct LoopsForInterpolation {
+ /* This field points to a loop data which is to be used for interpolation.
+	 * The idea is to avoid unnecessary allocations for regular faces, where
+	 * we can simply 
+	 */
+	const CustomData *loop_data;
+	/* Loops data calculated for ptex corners. There are always 4 elements
+	 * in this custom data, aligned the following way:
+	 *
+	 *   index 0 -> uv (0, 0)
+	 *   index 1 -> uv (0, 1)
+	 *   index 2 -> uv (1, 1)
+	 *   index 3 -> uv (1, 0)
+	 *
+	 * Is allocated for non-regular faces (triangles and n-gons).
+	 */
+	CustomData loop_data_storage;
+	bool loop_data_storage_allocated;
+	/* Infices within loop_data to interpolate for. The indices are aligned with
+	 * uv coordinates in a similar way as indices in loop_data_storage.
+	 */
+	int loop_indices[4];
+} LoopsForInterpolation;
+
+static void loop_interpolation_init(
+        const SubdivMeshContext *ctx,
+        LoopsForInterpolation *loop_interpolation,
+		const MPoly *coarse_poly)
+{
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	if (coarse_poly->totloop == 4) {
+		loop_interpolation->loop_data = &coarse_mesh->ldata;
+		loop_interpolation->loop_indices[0] = coarse_poly->loopstart + 0;
+		loop_interpolation->loop_indices[1] = coarse_poly->loopstart + 1;
+		loop_interpolation->loop_indices[2] = coarse_poly->loopstart + 2;
+		loop_interpolation->loop_indices[3] = coarse_poly->loopstart + 3;
+		loop_interpolation->loop_data_storage_allocated = false;
+	}
+	else {
+		loop_interpolation->loop_data = &loop_interpolation->loop_data_storage;
+		/* Allocate storage for loops corresponding to ptex corners. */
+		CustomData_copy(&ctx->coarse_mesh->ldata,
+		                &loop_interpolation->loop_data_storage,
+		                CD_MASK_EVERYTHING,
+		                CD_CALLOC,
+		                4);
+		/* Initialize indices. */
+		loop_interpolation->loop_indices[0] = 0;
+		loop_interpolation->loop_indices[1] = 1;
+		loop_interpolation->loop_indices[2] = 2;
+		loop_interpolation->loop_indices[3] = 3;
+		loop_interpolation->loop_data_storage_allocated = true;
+		/* Interpolate center of poly right away, it stays unchanged for all
+		 * ptex faces.
+		 */
+		const float weight = 1.0f / (float)coarse_poly->totloop;
+		float *weights = BLI_array_alloca(weights, coarse_poly->totloop);
+		int *indices = BLI_array_alloca(indices, coarse_poly->totloop);
+		for (int i = 0; i < coarse_poly->totloop; ++i) {
+			weights[i] = weight;
+			indices[i] = coarse_poly->loopstart + i;
+		}
+		CustomData_interp(&coarse_mesh->ldata,
+		                  &loop_interpolation->loop_data_storage,
+		                  indices,
+		                  weights, NULL,
+		                  coarse_poly->totloop,
+		                  2);
+	}
+}
+
+static void loop_interpolation_from_ptex(
+        const SubdivMeshContext *ctx,
+        LoopsForInterpolation *loop_interpolation,
+        const MPoly *coarse_poly,
+        const int ptex_face_index)
+{
+	if (coarse_poly->totloop == 4) {
+		/* Nothing to do, all indices and data is already assigned. */
+	} else {
+		const CustomData *loop_data = &ctx->coarse_mesh->ldata;
+		const Mesh *coarse_mesh = ctx->coarse_mesh;
+		const MLoop *coarse_mloop = coarse_mesh->mloop;
+		LoopsOfPtex loops_of_ptex;
+		loops_of_ptex_get(ctx, &loops_of_ptex, coarse_poly, ptex_face_index);
+		/* Ptex face corner corresponds to a poly loop with same index. */
+		CustomData_copy_data(loop_data,
+		                     &loop_interpolation->loop_data_storage,
+		                     coarse_poly->loopstart + ptex_face_index,
+		                     0,
+		                     1);
+		/* Interpolate remaining ptex face corners, which hits loops
+		 * middle points.
+		 *
+		 * TODO(sergey): Re-use one of interpolation results from previous
+		 * iteration.
+		 */
+		const float weights[2] = {0.5f, 0.5f};
+		const int first_indices[2] = {
+		        loops_of_ptex.first_loop - coarse_mloop,
+		        (loops_of_ptex.first_loop + 1 - coarse_mloop) %
+		                coarse_poly->totloop};
+		const int last_indices[2] = {
+		        loops_of_ptex.last_loop - coarse_mloop,
+		        loops_of_ptex.first_loop - coarse_mloop};
+		CustomData_interp(loop_data,
+		                  &loop_interpolation->loop_data_storage,
+		                  first_indices,
+		                  weights, NULL,
+		                  2,
+		                  1);
+		CustomData_interp(loop_data,
+		                  &loop_interpolation->loop_data_storage,
+		                  last_indices,
+		                  weights, NULL,
+		                  2,
+		                  3);
+	}
+}
+
+static void loop_interpolation_end(LoopsForInterpolation *loop_interpolation)
+{
+	if (loop_interpolation->loop_data_storage_allocated) {
+		CustomData_free(&loop_interpolation->loop_data_storage, 4);
+	}
+}
+
+static void subdiv_copy_vertex_data(
+        const SubdivMeshContext *ctx,
+        MVert *subdiv_vertex,
+        const VerticesForInterpolation *vertex_interpolation,
+        const float u, const float v)
+{
+	const int subdiv_vertex_index = subdiv_vertex - ctx->subdiv_mesh->mvert;
+	const float weights[4] = {(1.0f - u) * (1.0f - v),
+	                          u * (1.0f - v),
+	                          u * v,
+	                          (1.0f - u) * v};
+	CustomData_interp(vertex_interpolation->vertex_data,
+	                  &ctx->subdiv_mesh->vdata,
+	                  vertex_interpolation->vertex_indices,
+	                  weights, NULL,
+	                  4,
+	                  subdiv_vertex_index);
+	/* TODO(sergey): Set ORIGINDEX. */
+}
+
+static void subdiv_evaluate_vertices(SubdivMeshContext *ctx,
+                                     const int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MVert *subdiv_vertert = subdiv_mesh->mvert;
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	/* Actual evaluation. */
+	VerticesForInterpolation vertex_interpolation;
+	vertex_interpolation_init(ctx, &vertex_interpolation, coarse_poly);
+	MVert *subdiv_vert = &subdiv_vertert[ptex_face_index * resolution2];
+	for (int ptex_of_poly_index = 0;
+	     ptex_of_poly_index < num_poly_ptex_faces;
+	     ptex_of_poly_index++)
+	{
+		vertex_interpolation_from_ptex(ctx,
+		                               &vertex_interpolation,
+		                               coarse_poly,
+		                               ptex_of_poly_index);
+		const int current_ptex_face_index =
+		        ptex_face_index + ptex_of_poly_index;
+		BKE_subdiv_eval_limit_patch_resolution_point_and_short_normal(
+		        subdiv,
+		        current_ptex_face_index,
+		        resolution,
+		        subdiv_vert, offsetof(MVert, co), sizeof(MVert),
+		        subdiv_vert, offsetof(MVert, no), sizeof(MVert));
+		for (int y = 0; y < resolution; y++) {
+			const float v = y * inv_resolution_1;
+			for (int x = 0; x < resolution; x++, subdiv_vert++) {
+				const float u = x * inv_resolution_1;
+				subdiv_copy_vertex_data(ctx,
+				                        subdiv_vert,
+				                        &vertex_interpolation,
+				                        u, v);
+			}
+		}
+	}
+	vertex_interpolation_end(&vertex_interpolation);
+}
+
+static void subdiv_copy_edge_data(SubdivMeshContext *ctx,
+	                              MEdge *subdiv_edge,
+                                  const MEdge *coarse_edge)
+{
+	if (coarse_edge == NULL) {
+		subdiv_edge->crease = 0;
+		subdiv_edge->bweight = 0;
+		subdiv_edge->flag = 0;
+		return;
+	}
+	const int coarse_edge_index = coarse_edge - ctx->coarse_mesh->medge;
+	const int subdiv_edge_index = subdiv_edge - ctx->subdiv_mesh->medge;
+	CustomData_copy_data(&ctx->coarse_mesh->edata,
+	                     &ctx->subdiv_mesh->edata,
+	                     coarse_edge_index,
+	                     subdiv_edge_index,
+	                     1);
+	if (ctx->edge_origindex != NULL) {
+		ctx->edge_origindex[subdiv_edge_index] = coarse_edge_index;
+	}
+}
+
+static MEdge *subdiv_create_edges_row(SubdivMeshContext *ctx,
+                                      MEdge *subdiv_edge,
+                                      const MEdge *coarse_edge,
+                                      const int start_vertex_index,
+                                      const int resolution)
+{
+	int vertex_index = start_vertex_index;
+	for (int edge_index = 0;
+	     edge_index < resolution - 1;
+	     edge_index++, subdiv_edge++)
+	{
+		subdiv_copy_edge_data(ctx, subdiv_edge, coarse_edge);
+		subdiv_edge->v1 = vertex_index;
+		subdiv_edge->v2 = vertex_index + 1;
+		vertex_index += 1;
+	}
+	return subdiv_edge;
+}
+
+static MEdge *subdiv_create_edges_column(SubdivMeshContext *ctx,
+                                         MEdge *subdiv_edge,
+                                         const MEdge *coarse_start_edge,
+                                         const MEdge *coarse_end_edge,
+                                         const int start_vertex_index,
+                                         const int resolution)
+{
+	int vertex_index = start_vertex_index;
+	for (int edge_index = 0;
+	     edge_index < resolution;
+	     edge_index++, subdiv_edge++)
+	{
+		const MEdge *coarse_edge = NULL;
+		if (edge_index == 0) {
+			coarse_edge = coarse_start_edge;
+		}
+		else if (edge_index == resolution - 1) {
+			coarse_edge = coarse_end_edge;
+		}
+		subdiv_copy_edge_data(ctx, subdiv_edge, coarse_edge);
+		subdiv_edge->v1 = vertex_index;
+		subdiv_edge->v2 = vertex_index + resolution;
+		vertex_index += 1;
+	}
+	return subdiv_edge;
+}
+
+static void subdiv_create_edges(SubdivMeshContext *ctx, int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	const int num_edges_per_ptex = num_edges_per_ptex_get(resolution);
+	const int start_edge_index = ptex_face_index * num_edges_per_ptex;
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MEdge *subdiv_medge = subdiv_mesh->medge;
+	MEdge *subdiv_edge = &subdiv_medge[start_edge_index];
+	const int start_poly_vertex_index = ptex_face_index * resolution2;
+	/* Consider a subdivision of base face at level 1:
+	 *
+	 *  y
+	 *  ^
+	 *  |   (6) ---- (7) ---- (8)
+	 *  |    |        |        |
+	 *  |   (3) ---- (4) ---- (5)
+	 *  |    |        |        |
+	 *  |   (0) ---- (1) ---- (2)
+	 *  o---------------------------> x
+	 *
+	 * This is illustrate which parts of geometry is created by code below.
+	 */
+	for (int i = 0; i < num_poly_ptex_faces; i++) {
+		const int start_ptex_face_vertex_index =
+		        start_poly_vertex_index + i * resolution2;
+		EdgesOfPtex edges_of_ptex;
+		edges_of_ptex_get(ctx, &edges_of_ptex, coarse_poly, i);
+		/* Create bottom row of edges (0-1, 1-2). */
+		subdiv_edge = subdiv_create_edges_row(ctx,
+		                                      subdiv_edge,
+		                                      edges_of_ptex.first_edge,
+	                                          start_ptex_face_vertex_index,
+		                                      resolution);
+		/* Create remaining edges. */
+		for (int row = 0; row < resolution - 1; row++) {
+			const int start_row_vertex_index =
+			    start_ptex_face_vertex_index + row * resolution;
+			/* Create vertical columns.
+			 *
+			 * At first iteration it will be edges (0-3. 1-4, 2-5), then it
+			 * will be (3-6, 4-7, 5-8) and so on.
+			 */
+			subdiv_edge = subdiv_create_edges_column(
+			        ctx,
+			        subdiv_edge,
+			        edges_of_ptex.last_edge,
+			        edges_of_ptex.second_edge,
+			        start_row_vertex_index,
+			        resolution);
+			/* Create horizontal edge row.
+			 *
+			 * At first iteration it will be edges (3-4, 4-5), then it will be
+			 * (6-7, 7-8) and so on.
+			 */
+			subdiv_edge = subdiv_create_edges_row(
+			        ctx,
+			        subdiv_edge,
+			        (row == resolution - 2) ? edges_of_ptex.third_edge
+			                                : NULL,
+			        start_row_vertex_index + resolution,
+			        resolution);
+		}
+	}
+}
+
+static void subdiv_copy_loop_data(
+        const SubdivMeshContext *ctx,
+        MLoop *subdiv_loop,
+        const LoopsForInterpolation *loop_interpolation,
+        const float u, const float v)
+{
+	const int subdiv_loop_index = subdiv_loop - ctx->subdiv_mesh->mloop;
+	const float weights[4] = {(1.0f - u) * (1.0f - v),
+	                          u * (1.0f - v),
+	                          u * v,
+	                          (1.0f - u) * v};
+	CustomData_interp(loop_interpolation->loop_data,
+	                  &ctx->subdiv_mesh->ldata,
+	                  loop_interpolation->loop_indices,
+	                  weights, NULL,
+	                  4,
+	                  subdiv_loop_index);
+	/* TODO(sergey): Set ORIGINDEX. */
+}
+
+static void subdiv_eval_uv_layer(SubdivMeshContext *ctx,
+                                 MLoop *subdiv_loop,
+                                 const int ptex_face_index,
+                                 const float u, const float v,
+                                 const float inv_resolution_1)
+{
+	if (ctx->num_uv_layers == 0) {
+		return;
+	}
+	Subdiv *subdiv = ctx->subdiv;
+	const int mloop_index = subdiv_loop - ctx->subdiv_mesh->mloop;
+	const float du = inv_resolution_1;
+	const float dv = inv_resolution_1;
+	for (int layer_index = 0; layer_index < ctx->num_uv_layers; layer_index++) {
+		MLoopUV *subdiv_loopuv = &ctx->uv_layers[layer_index][mloop_index];
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u, v,
+		                             subdiv_loopuv[0].uv);
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u + du, v,
+		                             subdiv_loopuv[1].uv);
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u + du, v + dv,
+		                             subdiv_loopuv[2].uv);
+		BKE_subdiv_eval_face_varying(subdiv,
+		                             ptex_face_index,
+		                             u, v + dv,
+		                             subdiv_loopuv[3].uv);
+		/* TODO(sergey): Currently evaluator only has single UV layer, so can
+		 * not evaluate more than that. Need to be solved.
+		 */
+		break;
+	}
+}
+
+static void subdiv_create_loops(SubdivMeshContext *ctx, int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const float inv_resolution_1 = 1.0f / (float)(resolution - 1);
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	const int num_edges_per_ptex = num_edges_per_ptex_get(resolution);
+	const int start_edge_index = ptex_face_index * num_edges_per_ptex;
+	const int num_polys_per_ptex = num_polys_per_ptex_get(resolution);
+	const int start_poly_index = ptex_face_index * num_polys_per_ptex;
+	const int start_loop_index = 4 * start_poly_index;
+	const int start_vert_index = ptex_face_index * resolution2;
+	const float du = inv_resolution_1;
+	const float dv = inv_resolution_1;
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MLoop *subdiv_loopoop = subdiv_mesh->mloop;
+	MLoop *subdiv_loop = &subdiv_loopoop[start_loop_index];
+	LoopsForInterpolation loop_interpolation;
+	loop_interpolation_init(ctx, &loop_interpolation, coarse_poly);
+	for (int ptex_of_poly_index = 0;
+	     ptex_of_poly_index < num_poly_ptex_faces;
+	     ptex_of_poly_index++)
+	{
+		loop_interpolation_from_ptex(ctx,
+		                             &loop_interpolation,
+		                             coarse_poly,
+		                             ptex_of_poly_index);
+		const int current_ptex_face_index =
+		        ptex_face_index + ptex_of_poly_index;
+		for (int y = 0; y < resolution - 1; y++) {
+			const float v = y * inv_resolution_1;
+			for (int x = 0; x < resolution - 1; x++, subdiv_loop += 4) {
+				const float u = x * inv_resolution_1;
+				/* Vertex indicies ordered counter-clockwise. */
+				const int v0 = start_vert_index +
+				               (ptex_of_poly_index * resolution2) +
+				               (y * resolution + x);
+				const int v1 = v0 + 1;
+				const int v2 = v0 + resolution + 1;
+				const int v3 = v0 + resolution;
+				/* Edge indicies ordered counter-clockwise. */
+				const int e0 = start_edge_index +
+				               (ptex_of_poly_index * num_edges_per_ptex) +
+				               (y * (2 * resolution - 1) + x);
+				const int e1 = e0 + resolution;
+				const int e2 = e0 + (2 * resolution - 1);
+				const int e3 = e0 + resolution - 1;
+				/* Initialize 4 loops of corresponding hi-poly poly. */
+				/* TODO(sergey): For ptex boundaries we should use loops from
+				 * coarse mesh.
+				 */
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[0],
+				                      &loop_interpolation,
+				                      u, v);
+				subdiv_loop[0].v = v0;
+				subdiv_loop[0].e = e0;
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[1],
+				                      &loop_interpolation,
+				                      u + du, v);
+				subdiv_loop[1].v = v1;
+				subdiv_loop[1].e = e1;
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[2],
+				                      &loop_interpolation,
+				                      u + du, v + dv);
+				subdiv_loop[2].v = v2;
+				subdiv_loop[2].e = e2;
+				subdiv_copy_loop_data(ctx,
+				                      &subdiv_loop[3],
+				                      &loop_interpolation,
+				                      u, v + dv);
+				subdiv_loop[3].v = v3;
+				subdiv_loop[3].e = e3;
+				/* Interpolate UV layers using OpenSubdiv. */
+				subdiv_eval_uv_layer(ctx,
+				                     subdiv_loop,
+				                     current_ptex_face_index,
+				                     u, v,
+				                     inv_resolution_1);
+			}
+		}
+	}
+	loop_interpolation_end(&loop_interpolation);
+}
+
+static void subdiv_copy_poly_data(const SubdivMeshContext *ctx,
+                                  MPoly *subdiv_poly,
+                                  const MPoly *coarse_poly)
+{
+	const int coarse_poly_index = coarse_poly - ctx->coarse_mesh->mpoly;
+	const int subdiv_poly_index = subdiv_poly - ctx->subdiv_mesh->mpoly;
+	CustomData_copy_data(&ctx->coarse_mesh->pdata,
+	                     &ctx->subdiv_mesh->pdata,
+	                     coarse_poly_index,
+	                     subdiv_poly_index,
+	                     1);
+	if (ctx->poly_origindex != NULL) {
+		ctx->poly_origindex[subdiv_poly_index] = coarse_poly_index;
+	}
+}
+
+static void subdiv_create_polys(SubdivMeshContext *ctx, int poly_index)
+{
+	Subdiv *subdiv = ctx->subdiv;
+	const int resolution = ctx->settings->resolution;
+	const int ptex_face_index = subdiv->face_ptex_offset[poly_index];
+	const int num_polys_per_ptex = num_polys_per_ptex_get(resolution);
+	const int num_loops_per_ptex = 4 * num_polys_per_ptex;
+	const int start_poly_index = ptex_face_index * num_polys_per_ptex;
+	const int start_loop_index = 4 * start_poly_index;
+	/* Base/coarse mesh information. */
+	const Mesh *coarse_mesh = ctx->coarse_mesh;
+	const MPoly *coarse_polyoly = coarse_mesh->mpoly;
+	const MPoly *coarse_poly = &coarse_polyoly[poly_index];
+	const int num_poly_ptex_faces = mpoly_ptex_faces_count_get(coarse_poly);
+	/* Hi-poly subdivided mesh. */
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	MPoly *subdiv_mpoly = subdiv_mesh->mpoly;
+	MPoly *subdiv_mp = &subdiv_mpoly[start_poly_index];
+	for (int ptex_of_poly_index = 0;
+	     ptex_of_poly_index < num_poly_ptex_faces;
+	     ptex_of_poly_index++)
+	{
+		for (int subdiv_poly_index = 0;
+		     subdiv_poly_index < num_polys_per_ptex;
+		     subdiv_poly_index++, subdiv_mp++)
+		{
+			subdiv_copy_poly_data(ctx, subdiv_mp, coarse_poly);
+			subdiv_mp->loopstart = start_loop_index +
+			                       (ptex_of_poly_index * num_loops_per_ptex) +
+			                       (subdiv_poly_index * 4);
+			subdiv_mp->totloop = 4;
+		}
+	}
+}
+
+static void subdiv_eval_task(
+        void *__restrict userdata,
+        const int poly_index,
+        const ParallelRangeTLS *__restrict UNUSED(tls))
+{
+	SubdivMeshContext *data = userdata;
+	/* Evaluate hi-poly vertex coordinates and normals. */
+	subdiv_evaluate_vertices(data, poly_index);
+	/* Create mesh geometry for the given base poly index. */
+	subdiv_create_edges(data, poly_index);
+	subdiv_create_loops(data, poly_index);
+	subdiv_create_polys(data, poly_index);
+}
+
+static void cache_uv_layers(SubdivMeshContext *ctx)
+{
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	ctx->num_uv_layers =
+	        CustomData_number_of_layers(&subdiv_mesh->ldata, CD_MLOOPUV);
+	for (int layer_index = 0; layer_index < ctx->num_uv_layers; ++layer_index) {
+		ctx->uv_layers[layer_index] = CustomData_get_layer_n(
+		        &subdiv_mesh->ldata, CD_MLOOPUV, layer_index);
+	}
+}
+
+static void cache_custom_data_layers(SubdivMeshContext *ctx)
+{
+	Mesh *subdiv_mesh = ctx->subdiv_mesh;
+	/* Pointers to original indices layers. */
+	ctx->vert_origindex = CustomData_get_layer(
+	        &subdiv_mesh->vdata, CD_ORIGINDEX);
+	ctx->edge_origindex = CustomData_get_layer(
+	        &subdiv_mesh->edata, CD_ORIGINDEX);
+	ctx->loop_origindex = CustomData_get_layer(
+	        &subdiv_mesh->ldata, CD_ORIGINDEX);
+	ctx->poly_origindex = CustomData_get_layer(
+	        &subdiv_mesh->pdata, CD_ORIGINDEX);
+	/* UV layers interpolation. */
+	cache_uv_layers(ctx);
+}
+
+Mesh *BKE_subdiv_to_mesh(
+        Subdiv *subdiv,
+        const SubdivToMeshSettings *settings,
+        const Mesh *coarse_mesh)
+{
+	/* Make sure evaluator is up to date with possible new topology, and that
+	 * is is refined for the new positions of coarse vertices.
+	 */
+	BKE_subdiv_eval_update_from_mesh(subdiv, coarse_mesh);
+	const int resolution = settings->resolution;
+	const int resolution2 = resolution * resolution;
+	const int num_result_verts = subdiv->num_ptex_faces * resolution2;
+	const int num_result_edges =
+			subdiv->num_ptex_faces * num_edges_per_ptex_get(resolution);
+	const int num_result_polys =
+			subdiv->num_ptex_faces * num_polys_per_ptex_get(resolution);
+	const int num_result_loops = 4 * num_result_polys;
+	/* Create mesh and its arrays. */
+	Mesh *result = BKE_mesh_new_nomain_from_template(
+	        coarse_mesh,
+	        num_result_verts,
+	        num_result_edges,
+	        0,
+	        num_result_loops,
+	        num_result_polys);
+	/* Evaluate subdivisions of base faces in threads. */
+	SubdivMeshContext ctx;
+	ctx.coarse_mesh = coarse_mesh;
+	ctx.subdiv = subdiv;
+	ctx.subdiv_mesh = result;
+	ctx.settings = settings;
+	cache_custom_data_layers(&ctx);
+	/* Multi-threaded evaluation. */
+	ParallelRangeSettings parallel_range_settings;
+	BLI_parallel_range_settings_defaults(&parallel_range_settings);
+	BLI_task_parallel_range(0, coarse_mesh->totpoly,
+	                        &ctx,
+	                        subdiv_eval_task,
+	                        &parallel_range_settings);
+	return result;
+}