Solidify Modifier: support non-manifold input

This adds a new mode to solidify to support non-manifold geometry
with edges using 3 or more faces as input, resulting in a manifold mesh.

Since the differences between these methods don't translate well
into short terms, they're named "Simple" and "Complex" in the UI.

This also adds clamp with respect to angles
to the existing solidify modifier calculation.

5 files changed, 3426 insertions, 917 deletions

diff --git a/source/blender/modifiers/CMakeLists.txt b/source/blender/modifiers/CMakeLists.txt
index 1ae1f891e6f..602ef634e1b 100644
--- a/source/blender/modifiers/CMakeLists.txt
+++ b/source/blender/modifiers/CMakeLists.txt
@@ -85,6 +85,8 @@ set(SRC
   intern/MOD_smooth.c
   intern/MOD_softbody.c
   intern/MOD_solidify.c
+  intern/MOD_solidify_extrude.c
+  intern/MOD_solidify_nonmanifold.c
   intern/MOD_subsurf.c
   intern/MOD_surface.c
   intern/MOD_surfacedeform.c
@@ -104,6 +106,7 @@ set(SRC
   MOD_modifiertypes.h
   intern/MOD_fluidsim_util.h
   intern/MOD_meshcache_util.h
+  intern/MOD_solidify_util.h
   intern/MOD_util.h
   intern/MOD_weightvg_util.h
 )
diff --git a/source/blender/modifiers/intern/MOD_solidify.c b/source/blender/modifiers/intern/MOD_solidify.c
index 292e659fe03..8ea0a602b65 100644
--- a/source/blender/modifiers/intern/MOD_solidify.c
+++ b/source/blender/modifiers/intern/MOD_solidify.c
@@ -23,142 +23,25 @@
 
 #include "BLI_utildefines.h"
 
-#include "BLI_bitmap.h"
-#include "BLI_math.h"
-#include "BLI_utildefines_stack.h"
-
 #include "DNA_mesh_types.h"
-#include "DNA_meshdata_types.h"
-
-#include "MEM_guardedalloc.h"
 
-#include "BKE_mesh.h"
 #include "BKE_particle.h"
-#include "BKE_deform.h"
 
 #include "MOD_modifiertypes.h"
-#include "MOD_util.h"
+
+#include "MOD_solidify_util.h"
 
 #ifdef __GNUC__
 #  pragma GCC diagnostic error "-Wsign-conversion"
 #endif
 
-/* skip shell thickness for non-manifold edges, see [#35710] */
-#define USE_NONMANIFOLD_WORKAROUND
-
-/* *** derived mesh high quality normal calculation function  *** */
-/* could be exposed for other functions to use */
-
-typedef struct EdgeFaceRef {
-  int p1; /* init as -1 */
-  int p2;
-} EdgeFaceRef;
-
-BLI_INLINE bool edgeref_is_init(const EdgeFaceRef *edge_ref)
+static bool dependsOnNormals(ModifierData *md)
 {
-  return !((edge_ref->p1 == 0) && (edge_ref->p2 == 0));
-}
-
-/**
- * \param dm: Mesh to calculate normals for.
- * \param face_nors: Precalculated face normals.
- * \param r_vert_nors: Return vert normals.
- */
-static void mesh_calc_hq_normal(Mesh *mesh, float (*poly_nors)[3], float (*r_vert_nors)[3])
-{
-  int i, numVerts, numEdges, numPolys;
-  MPoly *mpoly, *mp;
-  MLoop *mloop, *ml;
-  MEdge *medge, *ed;
-  MVert *mvert, *mv;
-
-  numVerts = mesh->totvert;
-  numEdges = mesh->totedge;
-  numPolys = mesh->totpoly;
-  mpoly = mesh->mpoly;
-  medge = mesh->medge;
-  mvert = mesh->mvert;
-  mloop = mesh->mloop;
-
-  /* we don't want to overwrite any referenced layers */
-
-  /* Doesn't work here! */
-#if 0
-  mv = CustomData_duplicate_referenced_layer(&dm->vertData, CD_MVERT, numVerts);
-  cddm->mvert = mv;
-#endif
-
-  mv = mvert;
-  mp = mpoly;
-
-  {
-    EdgeFaceRef *edge_ref_array = MEM_calloc_arrayN(
-        (size_t)numEdges, sizeof(EdgeFaceRef), "Edge Connectivity");
-    EdgeFaceRef *edge_ref;
-    float edge_normal[3];
-
-    /* Add an edge reference if it's not there, pointing back to the face index. */
-    for (i = 0; i < numPolys; i++, mp++) {
-      int j;
-
-      ml = mloop + mp->loopstart;
-
-      for (j = 0; j < mp->totloop; j++, ml++) {
-        /* --- add edge ref to face --- */
-        edge_ref = &edge_ref_array[ml->e];
-        if (!edgeref_is_init(edge_ref)) {
-          edge_ref->p1 = i;
-          edge_ref->p2 = -1;
-        }
-        else if ((edge_ref->p1 != -1) && (edge_ref->p2 == -1)) {
-          edge_ref->p2 = i;
-        }
-        else {
-          /* 3+ faces using an edge, we can't handle this usefully */
-          edge_ref->p1 = edge_ref->p2 = -1;
-#ifdef USE_NONMANIFOLD_WORKAROUND
-          medge[ml->e].flag |= ME_EDGE_TMP_TAG;
-#endif
-        }
-        /* --- done --- */
-      }
-    }
-
-    for (i = 0, ed = medge, edge_ref = edge_ref_array; i < numEdges; i++, ed++, edge_ref++) {
-      /* Get the edge vert indices, and edge value (the face indices that use it) */
-
-      if (edgeref_is_init(edge_ref) && (edge_ref->p1 != -1)) {
-        if (edge_ref->p2 != -1) {
-          /* We have 2 faces using this edge, calculate the edges normal
-           * using the angle between the 2 faces as a weighting */
-#if 0
-          add_v3_v3v3(edge_normal, face_nors[edge_ref->f1], face_nors[edge_ref->f2]);
-          normalize_v3_length(
-              edge_normal,
-              angle_normalized_v3v3(face_nors[edge_ref->f1], face_nors[edge_ref->f2]));
-#else
-          mid_v3_v3v3_angle_weighted(
-              edge_normal, poly_nors[edge_ref->p1], poly_nors[edge_ref->p2]);
-#endif
-        }
-        else {
-          /* only one face attached to that edge */
-          /* an edge without another attached- the weight on this is undefined */
-          copy_v3_v3(edge_normal, poly_nors[edge_ref->p1]);
-        }
-        add_v3_v3(r_vert_nors[ed->v1], edge_normal);
-        add_v3_v3(r_vert_nors[ed->v2], edge_normal);
-      }
-    }
-    MEM_freeN(edge_ref_array);
-  }
-
-  /* normalize vertex normals and assign */
-  for (i = 0; i < numVerts; i++, mv++) {
-    if (normalize_v3(r_vert_nors[i]) == 0.0f) {
-      normal_short_to_float_v3(r_vert_nors[i], mv->no);
-    }
-  }
+  const SolidifyModifierData *smd = (SolidifyModifierData *)md;
+  /* even when we calculate our own normals,
+   * the vertex normals are used as a fallback
+   * if manifold is enabled vertex normals are not used */
+  return smd->mode == MOD_SOLIDIFY_MODE_EXTRUDE;
 }
 
 static void initData(ModifierData *md)
@@ -167,6 +50,9 @@ static void initData(ModifierData *md)
   smd->offset = 0.01f;
   smd->offset_fac = -1.0f;
   smd->flag = MOD_SOLIDIFY_RIM;
+  smd->mode = MOD_SOLIDIFY_MODE_EXTRUDE;
+  smd->nonmanifold_offset_mode = MOD_SOLIDIFY_NONMANIFOLD_OFFSET_MODE_CONSTRAINTS;
+  smd->nonmanifold_boundary_mode = MOD_SOLIDIFY_NONMANIFOLD_BOUNDARY_MODE_NONE;
 }
 
 static void requiredDataMask(Object *UNUSED(ob),
@@ -181,803 +67,16 @@ static void requiredDataMask(Object *UNUSED(ob),
   }
 }
 
-/* specific function for solidify - define locally */
-BLI_INLINE void madd_v3v3short_fl(float r[3], const short a[3], const float f)
-{
-  r[0] += (float)a[0] * f;
-  r[1] += (float)a[1] * f;
-  r[2] += (float)a[2] * f;
-}
-
 static Mesh *applyModifier(ModifierData *md, const ModifierEvalContext *ctx, Mesh *mesh)
 {
-  Mesh *result;
   const SolidifyModifierData *smd = (SolidifyModifierData *)md;
-
-  MVert *mv, *mvert, *orig_mvert;
-  MEdge *ed, *medge, *orig_medge;
-  MLoop *ml, *mloop, *orig_mloop;
-  MPoly *mp, *mpoly, *orig_mpoly;
-  const uint numVerts = (uint)mesh->totvert;
-  const uint numEdges = (uint)mesh->totedge;
-  const uint numPolys = (uint)mesh->totpoly;
-  const uint numLoops = (uint)mesh->totloop;
-  uint newLoops = 0, newPolys = 0, newEdges = 0, newVerts = 0, rimVerts = 0;
-
-  /* only use material offsets if we have 2 or more materials  */
-  const short mat_nr_max = ctx->object->totcol > 1 ? ctx->object->totcol - 1 : 0;
-  const short mat_ofs = mat_nr_max ? smd->mat_ofs : 0;
-  const short mat_ofs_rim = mat_nr_max ? smd->mat_ofs_rim : 0;
-
-  /* use for edges */
-  /* over-alloc new_vert_arr, old_vert_arr */
-  uint *new_vert_arr = NULL;
-  STACK_DECLARE(new_vert_arr);
-
-  uint *new_edge_arr = NULL;
-  STACK_DECLARE(new_edge_arr);
-
-  uint *old_vert_arr = MEM_calloc_arrayN(
-      numVerts, sizeof(*old_vert_arr), "old_vert_arr in solidify");
-
-  uint *edge_users = NULL;
-  char *edge_order = NULL;
-
-  float(*vert_nors)[3] = NULL;
-  float(*poly_nors)[3] = NULL;
-
-  const bool need_poly_normals = (smd->flag & MOD_SOLIDIFY_NORMAL_CALC) ||
-                                 (smd->flag & MOD_SOLIDIFY_EVEN);
-
-  const float ofs_orig = -(((-smd->offset_fac + 1.0f) * 0.5f) * smd->offset);
-  const float ofs_new = smd->offset + ofs_orig;
-  const float offset_fac_vg = smd->offset_fac_vg;
-  const float offset_fac_vg_inv = 1.0f - smd->offset_fac_vg;
-  const bool do_flip = (smd->flag & MOD_SOLIDIFY_FLIP) != 0;
-  const bool do_clamp = (smd->offset_clamp != 0.0f);
-  const bool do_shell = ((smd->flag & MOD_SOLIDIFY_RIM) && (smd->flag & MOD_SOLIDIFY_NOSHELL)) ==
-                        0;
-
-  /* weights */
-  MDeformVert *dvert;
-  const bool defgrp_invert = (smd->flag & MOD_SOLIDIFY_VGROUP_INV) != 0;
-  int defgrp_index;
-
-  /* array size is doubled in case of using a shell */
-  const uint stride = do_shell ? 2 : 1;
-
-  MOD_get_vgroup(ctx->object, mesh, smd->defgrp_name, &dvert, &defgrp_index);
-
-  orig_mvert = mesh->mvert;
-  orig_medge = mesh->medge;
-  orig_mloop = mesh->mloop;
-  orig_mpoly = mesh->mpoly;
-
-  if (need_poly_normals) {
-    /* calculate only face normals */
-    poly_nors = MEM_malloc_arrayN(numPolys, sizeof(*poly_nors), __func__);
-    BKE_mesh_calc_normals_poly(orig_mvert,
-                               NULL,
-                               (int)numVerts,
-                               orig_mloop,
-                               orig_mpoly,
-                               (int)numLoops,
-                               (int)numPolys,
-                               poly_nors,
-                               true);
-  }
-
-  STACK_INIT(new_vert_arr, numVerts * 2);
-  STACK_INIT(new_edge_arr, numEdges * 2);
-
-  if (smd->flag & MOD_SOLIDIFY_RIM) {
-    BLI_bitmap *orig_mvert_tag = BLI_BITMAP_NEW(numVerts, __func__);
-    uint eidx;
-    uint i;
-
-#define INVALID_UNUSED ((uint)-1)
-#define INVALID_PAIR ((uint)-2)
-
-    new_vert_arr = MEM_malloc_arrayN(numVerts, 2 * sizeof(*new_vert_arr), __func__);
-    new_edge_arr = MEM_malloc_arrayN(((numEdges * 2) + numVerts), sizeof(*new_edge_arr), __func__);
-
-    edge_users = MEM_malloc_arrayN(numEdges, sizeof(*edge_users), "solid_mod edges");
-    edge_order = MEM_malloc_arrayN(numEdges, sizeof(*edge_order), "solid_mod eorder");
-
-    /* save doing 2 loops here... */
-#if 0
-    copy_vn_i(edge_users, numEdges, INVALID_UNUSED);
-#endif
-
-    for (eidx = 0, ed = orig_medge; eidx < numEdges; eidx++, ed++) {
-      edge_users[eidx] = INVALID_UNUSED;
-    }
-
-    for (i = 0, mp = orig_mpoly; i < numPolys; i++, mp++) {
-      MLoop *ml_prev;
-      int j;
-
-      ml = orig_mloop + mp->loopstart;
-      ml_prev = ml + (mp->totloop - 1);
-
-      for (j = 0; j < mp->totloop; j++, ml++) {
-        /* add edge user */
-        eidx = ml_prev->e;
-        if (edge_users[eidx] == INVALID_UNUSED) {
-          ed = orig_medge + eidx;
-          BLI_assert(ELEM(ml_prev->v, ed->v1, ed->v2) && ELEM(ml->v, ed->v1, ed->v2));
-          edge_users[eidx] = (ml_prev->v > ml->v) == (ed->v1 < ed->v2) ? i : (i + numPolys);
-          edge_order[eidx] = j;
-        }
-        else {
-          edge_users[eidx] = INVALID_PAIR;
-        }
-        ml_prev = ml;
-      }
-    }
-
-    for (eidx = 0, ed = orig_medge; eidx < numEdges; eidx++, ed++) {
-      if (!ELEM(edge_users[eidx], INVALID_UNUSED, INVALID_PAIR)) {
-        BLI_BITMAP_ENABLE(orig_mvert_tag, ed->v1);
-        BLI_BITMAP_ENABLE(orig_mvert_tag, ed->v2);
-        STACK_PUSH(new_edge_arr, eidx);
-        newPolys++;
-        newLoops += 4;
-      }
-    }
-
-    for (i = 0; i < numVerts; i++) {
-      if (BLI_BITMAP_TEST(orig_mvert_tag, i)) {
-        old_vert_arr[i] = STACK_SIZE(new_vert_arr);
-        STACK_PUSH(new_vert_arr, i);
-        rimVerts++;
-      }
-      else {
-        old_vert_arr[i] = INVALID_UNUSED;
-      }
-    }
-
-    MEM_freeN(orig_mvert_tag);
-  }
-
-  if (do_shell == false) {
-    /* only add rim vertices */
-    newVerts = rimVerts;
-    /* each extruded face needs an opposite edge */
-    newEdges = newPolys;
-  }
-  else {
-    /* (stride == 2) in this case, so no need to add newVerts/newEdges */
-    BLI_assert(newVerts == 0);
-    BLI_assert(newEdges == 0);
-  }
-
-  if (smd->flag & MOD_SOLIDIFY_NORMAL_CALC) {
-    vert_nors = MEM_calloc_arrayN(numVerts, 3 * sizeof(float), "mod_solid_vno_hq");
-    mesh_calc_hq_normal(mesh, poly_nors, vert_nors);
-  }
-
-  result = BKE_mesh_new_nomain_from_template(mesh,
-                                             (int)((numVerts * stride) + newVerts),
-                                             (int)((numEdges * stride) + newEdges + rimVerts),
-                                             0,
-                                             (int)((numLoops * stride) + newLoops),
-                                             (int)((numPolys * stride) + newPolys));
-
-  mpoly = result->mpoly;
-  mloop = result->mloop;
-  medge = result->medge;
-  mvert = result->mvert;
-
-  if (do_shell) {
-    CustomData_copy_data(&mesh->vdata, &result->vdata, 0, 0, (int)numVerts);
-    CustomData_copy_data(&mesh->vdata, &result->vdata, 0, (int)numVerts, (int)numVerts);
-
-    CustomData_copy_data(&mesh->edata, &result->edata, 0, 0, (int)numEdges);
-    CustomData_copy_data(&mesh->edata, &result->edata, 0, (int)numEdges, (int)numEdges);
-
-    CustomData_copy_data(&mesh->ldata, &result->ldata, 0, 0, (int)numLoops);
-    /* DO NOT copy here the 'copied' part of loop data, we want to reverse loops
-     * (so that winding of copied face get reversed, so that normals get reversed
-     * and point in expected direction...).
-     * If we also copy data here, then this data get overwritten
-     * (and allocated memory becomes memleak). */
-
-    CustomData_copy_data(&mesh->pdata, &result->pdata, 0, 0, (int)numPolys);
-    CustomData_copy_data(&mesh->pdata, &result->pdata, 0, (int)numPolys, (int)numPolys);
+  if (smd->mode == MOD_SOLIDIFY_MODE_EXTRUDE) {
+    return MOD_solidify_extrude_applyModifier(md, ctx, mesh);
   }
-  else {
-    int i, j;
-    CustomData_copy_data(&mesh->vdata, &result->vdata, 0, 0, (int)numVerts);
-    for (i = 0, j = (int)numVerts; i < numVerts; i++) {
-      if (old_vert_arr[i] != INVALID_UNUSED) {
-        CustomData_copy_data(&mesh->vdata, &result->vdata, i, j, 1);
-        j++;
-      }
-    }
-
-    CustomData_copy_data(&mesh->edata, &result->edata, 0, 0, (int)numEdges);
-
-    for (i = 0, j = (int)numEdges; i < numEdges; i++) {
-      if (!ELEM(edge_users[i], INVALID_UNUSED, INVALID_PAIR)) {
-        MEdge *ed_src, *ed_dst;
-        CustomData_copy_data(&mesh->edata, &result->edata, i, j, 1);
-
-        ed_src = &medge[i];
-        ed_dst = &medge[j];
-        ed_dst->v1 = old_vert_arr[ed_src->v1] + numVerts;
-        ed_dst->v2 = old_vert_arr[ed_src->v2] + numVerts;
-        j++;
-      }
-    }
-
-    /* will be created later */
-    CustomData_copy_data(&mesh->ldata, &result->ldata, 0, 0, (int)numLoops);
-    CustomData_copy_data(&mesh->pdata, &result->pdata, 0, 0, (int)numPolys);
+  else if (smd->mode == MOD_SOLIDIFY_MODE_NONMANIFOLD) {
+    return MOD_solidify_nonmanifold_applyModifier(md, ctx, mesh);
   }
-
-#undef INVALID_UNUSED
-#undef INVALID_PAIR
-
-  /* initializes: (i_end, do_shell_align, mv)  */
-#define INIT_VERT_ARRAY_OFFSETS(test) \
-  if (((ofs_new >= ofs_orig) == do_flip) == test) { \
-    i_end = numVerts; \
-    do_shell_align = true; \
-    mv = mvert; \
-  } \
-  else { \
-    if (do_shell) { \
-      i_end = numVerts; \
-      do_shell_align = true; \
-    } \
-    else { \
-      i_end = newVerts; \
-      do_shell_align = false; \
-    } \
-    mv = &mvert[numVerts]; \
-  } \
-  (void)0
-
-  /* flip normals */
-
-  if (do_shell) {
-    uint i;
-
-    mp = mpoly + numPolys;
-    for (i = 0; i < mesh->totpoly; i++, mp++) {
-      const int loop_end = mp->totloop - 1;
-      MLoop *ml2;
-      uint e;
-      int j;
-
-      /* reverses the loop direction (MLoop.v as well as custom-data)
-       * MLoop.e also needs to be corrected too, done in a separate loop below. */
-      ml2 = mloop + mp->loopstart + mesh->totloop;
-#if 0
-      for (j = 0; j < mp->totloop; j++) {
-        CustomData_copy_data(&mesh->ldata,
-                             &result->ldata,
-                             mp->loopstart + j,
-                             mp->loopstart + (loop_end - j) + mesh->totloop,
-                             1);
-      }
-#else
-      /* slightly more involved, keep the first vertex the same for the copy,
-       * ensures the diagonals in the new face match the original. */
-      j = 0;
-      for (int j_prev = loop_end; j < mp->totloop; j_prev = j++) {
-        CustomData_copy_data(&mesh->ldata,
-                             &result->ldata,
-                             mp->loopstart + j,
-                             mp->loopstart + (loop_end - j_prev) + mesh->totloop,
-                             1);
-      }
-#endif
-
-      if (mat_ofs) {
-        mp->mat_nr += mat_ofs;
-        CLAMP(mp->mat_nr, 0, mat_nr_max);
-      }
-
-      e = ml2[0].e;
-      for (j = 0; j < loop_end; j++) {
-        ml2[j].e = ml2[j + 1].e;
-      }
-      ml2[loop_end].e = e;
-
-      mp->loopstart += mesh->totloop;
-
-      for (j = 0; j < mp->totloop; j++) {
-        ml2[j].e += numEdges;
-        ml2[j].v += numVerts;
-      }
-    }
-
-    for (i = 0, ed = medge + numEdges; i < numEdges; i++, ed++) {
-      ed->v1 += numVerts;
-      ed->v2 += numVerts;
-    }
-  }
-
-  /* note, copied vertex layers don't have flipped normals yet. do this after applying offset */
-  if ((smd->flag & MOD_SOLIDIFY_EVEN) == 0) {
-    /* no even thickness, very simple */
-    float scalar_short;
-    float scalar_short_vgroup;
-
-    /* for clamping */
-    float *vert_lens = NULL;
-    const float offset = fabsf(smd->offset) * smd->offset_clamp;
-    const float offset_sq = offset * offset;
-
-    if (do_clamp) {
-      uint i;
-
-      vert_lens = MEM_malloc_arrayN(numVerts, sizeof(float), "vert_lens");
-      copy_vn_fl(vert_lens, (int)numVerts, FLT_MAX);
-      for (i = 0; i < numEdges; i++) {
-        const float ed_len_sq = len_squared_v3v3(mvert[medge[i].v1].co, mvert[medge[i].v2].co);
-        vert_lens[medge[i].v1] = min_ff(vert_lens[medge[i].v1], ed_len_sq);
-        vert_lens[medge[i].v2] = min_ff(vert_lens[medge[i].v2], ed_len_sq);
-      }
-    }
-
-    if (ofs_new != 0.0f) {
-      uint i_orig, i_end;
-      bool do_shell_align;
-
-      scalar_short = scalar_short_vgroup = ofs_new / 32767.0f;
-
-      INIT_VERT_ARRAY_OFFSETS(false);
-
-      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
-        const uint i = do_shell_align ? i_orig : new_vert_arr[i_orig];
-        if (dvert) {
-          MDeformVert *dv = &dvert[i];
-          if (defgrp_invert) {
-            scalar_short_vgroup = 1.0f - defvert_find_weight(dv, defgrp_index);
-          }
-          else {
-            scalar_short_vgroup = defvert_find_weight(dv, defgrp_index);
-          }
-          scalar_short_vgroup = (offset_fac_vg + (scalar_short_vgroup * offset_fac_vg_inv)) *
-                                scalar_short;
-        }
-        if (do_clamp) {
-          /* always reset becaise we may have set before */
-          if (dvert == NULL) {
-            scalar_short_vgroup = scalar_short;
-          }
-          if (vert_lens[i] < offset_sq) {
-            float scalar = sqrtf(vert_lens[i]) / offset;
-            scalar_short_vgroup *= scalar;
-          }
-        }
-        madd_v3v3short_fl(mv->co, mv->no, scalar_short_vgroup);
-      }
-    }
-
-    if (ofs_orig != 0.0f) {
-      uint i_orig, i_end;
-      bool do_shell_align;
-
-      scalar_short = scalar_short_vgroup = ofs_orig / 32767.0f;
-
-      /* as above but swapped */
-      INIT_VERT_ARRAY_OFFSETS(true);
-
-      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
-        const uint i = do_shell_align ? i_orig : new_vert_arr[i_orig];
-        if (dvert) {
-          MDeformVert *dv = &dvert[i];
-          if (defgrp_invert) {
-            scalar_short_vgroup = 1.0f - defvert_find_weight(dv, defgrp_index);
-          }
-          else {
-            scalar_short_vgroup = defvert_find_weight(dv, defgrp_index);
-          }
-          scalar_short_vgroup = (offset_fac_vg + (scalar_short_vgroup * offset_fac_vg_inv)) *
-                                scalar_short;
-        }
-        if (do_clamp) {
-          /* always reset becaise we may have set before */
-          if (dvert == NULL) {
-            scalar_short_vgroup = scalar_short;
-          }
-          if (vert_lens[i] < offset_sq) {
-            float scalar = sqrtf(vert_lens[i]) / offset;
-            scalar_short_vgroup *= scalar;
-          }
-        }
-        madd_v3v3short_fl(mv->co, mv->no, scalar_short_vgroup);
-      }
-    }
-
-    if (do_clamp) {
-      MEM_freeN(vert_lens);
-    }
-  }
-  else {
-#ifdef USE_NONMANIFOLD_WORKAROUND
-    const bool check_non_manifold = (smd->flag & MOD_SOLIDIFY_NORMAL_CALC) != 0;
-#endif
-    /* same as EM_solidify() in editmesh_lib.c */
-    float *vert_angles = MEM_calloc_arrayN(
-        numVerts, 2 * sizeof(float), "mod_solid_pair"); /* 2 in 1 */
-    float *vert_accum = vert_angles + numVerts;
-    uint vidx;
-    uint i;
-
-    if (vert_nors == NULL) {
-      vert_nors = MEM_malloc_arrayN(numVerts, 3 * sizeof(float), "mod_solid_vno");
-      for (i = 0, mv = mvert; i < numVerts; i++, mv++) {
-        normal_short_to_float_v3(vert_nors[i], mv->no);
-      }
-    }
-
-    for (i = 0, mp = mpoly; i < numPolys; i++, mp++) {
-      /* #BKE_mesh_calc_poly_angles logic is inlined here */
-      float nor_prev[3];
-      float nor_next[3];
-
-      int i_curr = mp->totloop - 1;
-      int i_next = 0;
-
-      ml = &mloop[mp->loopstart];
-
-      sub_v3_v3v3(nor_prev, mvert[ml[i_curr - 1].v].co, mvert[ml[i_curr].v].co);
-      normalize_v3(nor_prev);
-
-      while (i_next < mp->totloop) {
-        float angle;
-        sub_v3_v3v3(nor_next, mvert[ml[i_curr].v].co, mvert[ml[i_next].v].co);
-        normalize_v3(nor_next);
-        angle = angle_normalized_v3v3(nor_prev, nor_next);
-
-        /* --- not related to angle calc --- */
-        if (angle < FLT_EPSILON) {
-          angle = FLT_EPSILON;
-        }
-
-        vidx = ml[i_curr].v;
-        vert_accum[vidx] += angle;
-
-#ifdef USE_NONMANIFOLD_WORKAROUND
-        /* skip 3+ face user edges */
-        if ((check_non_manifold == false) ||
-            LIKELY(((orig_medge[ml[i_curr].e].flag & ME_EDGE_TMP_TAG) == 0) &&
-                   ((orig_medge[ml[i_next].e].flag & ME_EDGE_TMP_TAG) == 0))) {
-          vert_angles[vidx] += shell_v3v3_normalized_to_dist(vert_nors[vidx], poly_nors[i]) *
-                               angle;
-        }
-        else {
-          vert_angles[vidx] += angle;
-        }
-#else
-        vert_angles[vidx] += shell_v3v3_normalized_to_dist(vert_nors[vidx], poly_nors[i]) * angle;
-#endif
-        /* --- end non-angle-calc section --- */
-
-        /* step */
-        copy_v3_v3(nor_prev, nor_next);
-        i_curr = i_next;
-        i_next++;
-      }
-    }
-
-    /* vertex group support */
-    if (dvert) {
-      MDeformVert *dv = dvert;
-      float scalar;
-
-      if (defgrp_invert) {
-        for (i = 0; i < numVerts; i++, dv++) {
-          scalar = 1.0f - defvert_find_weight(dv, defgrp_index);
-          scalar = offset_fac_vg + (scalar * offset_fac_vg_inv);
-          vert_angles[i] *= scalar;
-        }
-      }
-      else {
-        for (i = 0; i < numVerts; i++, dv++) {
-          scalar = defvert_find_weight(dv, defgrp_index);
-          scalar = offset_fac_vg + (scalar * offset_fac_vg_inv);
-          vert_angles[i] *= scalar;
-        }
-      }
-    }
-
-    if (do_clamp) {
-      float *vert_lens_sq = MEM_malloc_arrayN(numVerts, sizeof(float), "vert_lens");
-      const float offset = fabsf(smd->offset) * smd->offset_clamp;
-      const float offset_sq = offset * offset;
-      copy_vn_fl(vert_lens_sq, (int)numVerts, FLT_MAX);
-      for (i = 0; i < numEdges; i++) {
-        const float ed_len = len_squared_v3v3(mvert[medge[i].v1].co, mvert[medge[i].v2].co);
-        vert_lens_sq[medge[i].v1] = min_ff(vert_lens_sq[medge[i].v1], ed_len);
-        vert_lens_sq[medge[i].v2] = min_ff(vert_lens_sq[medge[i].v2], ed_len);
-      }
-      for (i = 0; i < numVerts; i++) {
-        if (vert_lens_sq[i] < offset_sq) {
-          float scalar = sqrtf(vert_lens_sq[i]) / offset;
-          vert_angles[i] *= scalar;
-        }
-      }
-      MEM_freeN(vert_lens_sq);
-    }
-
-    if (ofs_new != 0.0f) {
-      uint i_orig, i_end;
-      bool do_shell_align;
-
-      INIT_VERT_ARRAY_OFFSETS(false);
-
-      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
-        const uint i_other = do_shell_align ? i_orig : new_vert_arr[i_orig];
-        if (vert_accum[i_other]) { /* zero if unselected */
-          madd_v3_v3fl(
-              mv->co, vert_nors[i_other], ofs_new * (vert_angles[i_other] / vert_accum[i_other]));
-        }
-      }
-    }
-
-    if (ofs_orig != 0.0f) {
-      uint i_orig, i_end;
-      bool do_shell_align;
-
-      /* same as above but swapped, intentional use of 'ofs_new' */
-      INIT_VERT_ARRAY_OFFSETS(true);
-
-      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
-        const uint i_other = do_shell_align ? i_orig : new_vert_arr[i_orig];
-        if (vert_accum[i_other]) { /* zero if unselected */
-          madd_v3_v3fl(
-              mv->co, vert_nors[i_other], ofs_orig * (vert_angles[i_other] / vert_accum[i_other]));
-        }
-      }
-    }
-
-    MEM_freeN(vert_angles);
-  }
-
-  if (vert_nors) {
-    MEM_freeN(vert_nors);
-  }
-
-  /* must recalculate normals with vgroups since they can displace unevenly [#26888] */
-  if ((mesh->runtime.cd_dirty_vert & CD_MASK_NORMAL) || (smd->flag & MOD_SOLIDIFY_RIM) || dvert) {
-    result->runtime.cd_dirty_vert |= CD_MASK_NORMAL;
-  }
-  else if (do_shell) {
-    uint i;
-    /* flip vertex normals for copied verts */
-    mv = mvert + numVerts;
-    for (i = 0; i < numVerts; i++, mv++) {
-      negate_v3_short(mv->no);
-    }
-  }
-
-  if (smd->flag & MOD_SOLIDIFY_RIM) {
-    uint i;
-
-    /* bugger, need to re-calculate the normals for the new edge faces.
-     * This could be done in many ways, but probably the quickest way
-     * is to calculate the average normals for side faces only.
-     * Then blend them with the normals of the edge verts.
-     *
-     * at the moment its easiest to allocate an entire array for every vertex,
-     * even though we only need edge verts - campbell
-     */
-
-#define SOLIDIFY_SIDE_NORMALS
-
-#ifdef SOLIDIFY_SIDE_NORMALS
-    /* Note that, due to the code setting cd_dirty_vert a few lines above,
-     * do_side_normals is always false. - Sybren */
-    const bool do_side_normals = !(result->runtime.cd_dirty_vert & CD_MASK_NORMAL);
-    /* annoying to allocate these since we only need the edge verts, */
-    float(*edge_vert_nos)[3] = do_side_normals ?
-                                   MEM_calloc_arrayN(numVerts, 3 * sizeof(float), __func__) :
-                                   NULL;
-    float nor[3];
-#endif
-    const uchar crease_rim = smd->crease_rim * 255.0f;
-    const uchar crease_outer = smd->crease_outer * 255.0f;
-    const uchar crease_inner = smd->crease_inner * 255.0f;
-
-    int *origindex_edge;
-    int *orig_ed;
-    uint j;
-
-    if (crease_rim || crease_outer || crease_inner) {
-      result->cd_flag |= ME_CDFLAG_EDGE_CREASE;
-    }
-
-    /* add faces & edges */
-    origindex_edge = CustomData_get_layer(&result->edata, CD_ORIGINDEX);
-    orig_ed = (origindex_edge) ? &origindex_edge[(numEdges * stride) + newEdges] : NULL;
-    ed = &medge[(numEdges * stride) + newEdges]; /* start after copied edges */
-    for (i = 0; i < rimVerts; i++, ed++) {
-      ed->v1 = new_vert_arr[i];
-      ed->v2 = (do_shell ? new_vert_arr[i] : i) + numVerts;
-      ed->flag |= ME_EDGEDRAW | ME_EDGERENDER;
-
-      if (orig_ed) {
-        *orig_ed = ORIGINDEX_NONE;
-        orig_ed++;
-      }
-
-      if (crease_rim) {
-        ed->crease = crease_rim;
-      }
-    }
-
-    /* faces */
-    mp = mpoly + (numPolys * stride);
-    ml = mloop + (numLoops * stride);
-    j = 0;
-    for (i = 0; i < newPolys; i++, mp++) {
-      uint eidx = new_edge_arr[i];
-      uint pidx = edge_users[eidx];
-      int k1, k2;
-      bool flip;
-
-      if (pidx >= numPolys) {
-        pidx -= numPolys;
-        flip = true;
-      }
-      else {
-        flip = false;
-      }
-
-      ed = medge + eidx;
-
-      /* copy most of the face settings */
-      CustomData_copy_data(
-          &mesh->pdata, &result->pdata, (int)pidx, (int)((numPolys * stride) + i), 1);
-      mp->loopstart = (int)(j + (numLoops * stride));
-      mp->flag = mpoly[pidx].flag;
-
-      /* notice we use 'mp->totloop' which is later overwritten,
-       * we could lookup the original face but there's no point since this is a copy
-       * and will have the same value, just take care when changing order of assignment */
-
-      /* prev loop */
-      k1 = mpoly[pidx].loopstart + (((edge_order[eidx] - 1) + mp->totloop) % mp->totloop);
-
-      k2 = mpoly[pidx].loopstart + (edge_order[eidx]);
-
-      mp->totloop = 4;
-
-      CustomData_copy_data(
-          &mesh->ldata, &result->ldata, k2, (int)((numLoops * stride) + j + 0), 1);
-      CustomData_copy_data(
-          &mesh->ldata, &result->ldata, k1, (int)((numLoops * stride) + j + 1), 1);
-      CustomData_copy_data(
-          &mesh->ldata, &result->ldata, k1, (int)((numLoops * stride) + j + 2), 1);
-      CustomData_copy_data(
-          &mesh->ldata, &result->ldata, k2, (int)((numLoops * stride) + j + 3), 1);
-
-      if (flip == false) {
-        ml[j].v = ed->v1;
-        ml[j++].e = eidx;
-
-        ml[j].v = ed->v2;
-        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v2] + newEdges;
-
-        ml[j].v = (do_shell ? ed->v2 : old_vert_arr[ed->v2]) + numVerts;
-        ml[j++].e = (do_shell ? eidx : i) + numEdges;
-
-        ml[j].v = (do_shell ? ed->v1 : old_vert_arr[ed->v1]) + numVerts;
-        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v1] + newEdges;
-      }
-      else {
-        ml[j].v = ed->v2;
-        ml[j++].e = eidx;
-
-        ml[j].v = ed->v1;
-        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v1] + newEdges;
-
-        ml[j].v = (do_shell ? ed->v1 : old_vert_arr[ed->v1]) + numVerts;
-        ml[j++].e = (do_shell ? eidx : i) + numEdges;
-
-        ml[j].v = (do_shell ? ed->v2 : old_vert_arr[ed->v2]) + numVerts;
-        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v2] + newEdges;
-      }
-
-      if (origindex_edge) {
-        origindex_edge[ml[j - 3].e] = ORIGINDEX_NONE;
-        origindex_edge[ml[j - 1].e] = ORIGINDEX_NONE;
-      }
-
-      /* use the next material index if option enabled */
-      if (mat_ofs_rim) {
-        mp->mat_nr += mat_ofs_rim;
-        CLAMP(mp->mat_nr, 0, mat_nr_max);
-      }
-      if (crease_outer) {
-        /* crease += crease_outer; without wrapping */
-        char *cr = &(ed->crease);
-        int tcr = *cr + crease_outer;
-        *cr = tcr > 255 ? 255 : tcr;
-      }
-
-      if (crease_inner) {
-        /* crease += crease_inner; without wrapping */
-        char *cr = &(medge[numEdges + (do_shell ? eidx : i)].crease);
-        int tcr = *cr + crease_inner;
-        *cr = tcr > 255 ? 255 : tcr;
-      }
-
-#ifdef SOLIDIFY_SIDE_NORMALS
-      if (do_side_normals) {
-        normal_quad_v3(nor,
-                       mvert[ml[j - 4].v].co,
-                       mvert[ml[j - 3].v].co,
-                       mvert[ml[j - 2].v].co,
-                       mvert[ml[j - 1].v].co);
-
-        add_v3_v3(edge_vert_nos[ed->v1], nor);
-        add_v3_v3(edge_vert_nos[ed->v2], nor);
-      }
-#endif
-    }
-
-#ifdef SOLIDIFY_SIDE_NORMALS
-    if (do_side_normals) {
-      const MEdge *ed_orig = medge;
-      ed = medge + (numEdges * stride);
-      for (i = 0; i < rimVerts; i++, ed++, ed_orig++) {
-        float nor_cpy[3];
-        short *nor_short;
-        int k;
-
-        /* note, only the first vertex (lower half of the index) is calculated */
-        BLI_assert(ed->v1 < numVerts);
-        normalize_v3_v3(nor_cpy, edge_vert_nos[ed_orig->v1]);
-
-        for (k = 0; k < 2; k++) { /* loop over both verts of the edge */
-          nor_short = mvert[*(&ed->v1 + k)].no;
-          normal_short_to_float_v3(nor, nor_short);
-          add_v3_v3(nor, nor_cpy);
-          normalize_v3(nor);
-          normal_float_to_short_v3(nor_short, nor);
-        }
-      }
-
-      MEM_freeN(edge_vert_nos);
-    }
-#endif
-
-    MEM_freeN(new_vert_arr);
-    MEM_freeN(new_edge_arr);
-
-    MEM_freeN(edge_users);
-    MEM_freeN(edge_order);
-  }
-
-  if (old_vert_arr) {
-    MEM_freeN(old_vert_arr);
-  }
-
-  if (poly_nors) {
-    MEM_freeN(poly_nors);
-  }
-
-  if (numPolys == 0 && numEdges != 0) {
-    modifier_setError(md, "Faces needed for useful output");
-  }
-
-  return result;
-}
-
-#undef SOLIDIFY_SIDE_NORMALS
-
-static bool dependsOnNormals(ModifierData *UNUSED(md))
-{
-  /* even when we calculate our own normals,
-   * the vertex normals are used as a fallback */
-  return true;
+  return mesh;
 }
 
 ModifierTypeInfo modifierType_Solidify = {
diff --git a/source/blender/modifiers/intern/MOD_solidify_extrude.c b/source/blender/modifiers/intern/MOD_solidify_extrude.c
new file mode 100644
index 00000000000..8c5a2551c0b
--- /dev/null
+++ b/source/blender/modifiers/intern/MOD_solidify_extrude.c
@@ -0,0 +1,1105 @@
+/*
+ * 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.
+ */
+
+/** \file
+ * \ingroup modifiers
+ */
+
+#include "BLI_utildefines.h"
+
+#include "BLI_bitmap.h"
+#include "BLI_math.h"
+#include "BLI_utildefines_stack.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_mesh.h"
+#include "BKE_particle.h"
+#include "BKE_deform.h"
+
+#include "MOD_modifiertypes.h"
+#include "MOD_util.h"
+#include "MOD_solidify_util.h" /* own include */
+
+#ifdef __GNUC__
+#  pragma GCC diagnostic error "-Wsign-conversion"
+#endif
+
+/* -------------------------------------------------------------------- */
+/** \name Local Utilities
+ * \{ */
+
+/* specific function for solidify - define locally */
+BLI_INLINE void madd_v3v3short_fl(float r[3], const short a[3], const float f)
+{
+  r[0] += (float)a[0] * f;
+  r[1] += (float)a[1] * f;
+  r[2] += (float)a[2] * f;
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name High Quality Normal Calculation Function
+ * \{ */
+
+/* skip shell thickness for non-manifold edges, see [#35710] */
+#define USE_NONMANIFOLD_WORKAROUND
+
+/* *** derived mesh high quality normal calculation function  *** */
+/* could be exposed for other functions to use */
+
+typedef struct EdgeFaceRef {
+  int p1; /* init as -1 */
+  int p2;
+} EdgeFaceRef;
+
+BLI_INLINE bool edgeref_is_init(const EdgeFaceRef *edge_ref)
+{
+  return !((edge_ref->p1 == 0) && (edge_ref->p2 == 0));
+}
+
+/**
+ * \param dm: Mesh to calculate normals for.
+ * \param poly_nors: Precalculated face normals.
+ * \param r_vert_nors: Return vert normals.
+ */
+static void mesh_calc_hq_normal(Mesh *mesh, float (*poly_nors)[3], float (*r_vert_nors)[3])
+{
+  int i, numVerts, numEdges, numPolys;
+  MPoly *mpoly, *mp;
+  MLoop *mloop, *ml;
+  MEdge *medge, *ed;
+  MVert *mvert, *mv;
+
+  numVerts = mesh->totvert;
+  numEdges = mesh->totedge;
+  numPolys = mesh->totpoly;
+  mpoly = mesh->mpoly;
+  medge = mesh->medge;
+  mvert = mesh->mvert;
+  mloop = mesh->mloop;
+
+  /* we don't want to overwrite any referenced layers */
+
+  /* Doesn't work here! */
+#if 0
+  mv = CustomData_duplicate_referenced_layer(&dm->vertData, CD_MVERT, numVerts);
+  cddm->mvert = mv;
+#endif
+
+  mv = mvert;
+  mp = mpoly;
+
+  {
+    EdgeFaceRef *edge_ref_array = MEM_calloc_arrayN(
+        (size_t)numEdges, sizeof(EdgeFaceRef), "Edge Connectivity");
+    EdgeFaceRef *edge_ref;
+    float edge_normal[3];
+
+    /* Add an edge reference if it's not there, pointing back to the face index. */
+    for (i = 0; i < numPolys; i++, mp++) {
+      int j;
+
+      ml = mloop + mp->loopstart;
+
+      for (j = 0; j < mp->totloop; j++, ml++) {
+        /* --- add edge ref to face --- */
+        edge_ref = &edge_ref_array[ml->e];
+        if (!edgeref_is_init(edge_ref)) {
+          edge_ref->p1 = i;
+          edge_ref->p2 = -1;
+        }
+        else if ((edge_ref->p1 != -1) && (edge_ref->p2 == -1)) {
+          edge_ref->p2 = i;
+        }
+        else {
+          /* 3+ faces using an edge, we can't handle this usefully */
+          edge_ref->p1 = edge_ref->p2 = -1;
+#ifdef USE_NONMANIFOLD_WORKAROUND
+          medge[ml->e].flag |= ME_EDGE_TMP_TAG;
+#endif
+        }
+        /* --- done --- */
+      }
+    }
+
+    for (i = 0, ed = medge, edge_ref = edge_ref_array; i < numEdges; i++, ed++, edge_ref++) {
+      /* Get the edge vert indices, and edge value (the face indices that use it) */
+
+      if (edgeref_is_init(edge_ref) && (edge_ref->p1 != -1)) {
+        if (edge_ref->p2 != -1) {
+          /* We have 2 faces using this edge, calculate the edges normal
+           * using the angle between the 2 faces as a weighting */
+#if 0
+          add_v3_v3v3(edge_normal, face_nors[edge_ref->f1], face_nors[edge_ref->f2]);
+          normalize_v3_length(
+              edge_normal,
+              angle_normalized_v3v3(face_nors[edge_ref->f1], face_nors[edge_ref->f2]));
+#else
+          mid_v3_v3v3_angle_weighted(
+              edge_normal, poly_nors[edge_ref->p1], poly_nors[edge_ref->p2]);
+#endif
+        }
+        else {
+          /* only one face attached to that edge */
+          /* an edge without another attached- the weight on this is undefined */
+          copy_v3_v3(edge_normal, poly_nors[edge_ref->p1]);
+        }
+        add_v3_v3(r_vert_nors[ed->v1], edge_normal);
+        add_v3_v3(r_vert_nors[ed->v2], edge_normal);
+      }
+    }
+    MEM_freeN(edge_ref_array);
+  }
+
+  /* normalize vertex normals and assign */
+  for (i = 0; i < numVerts; i++, mv++) {
+    if (normalize_v3(r_vert_nors[i]) == 0.0f) {
+      normal_short_to_float_v3(r_vert_nors[i], mv->no);
+    }
+  }
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Main Solidify Function
+ * \{ */
+
+Mesh *MOD_solidify_extrude_applyModifier(ModifierData *md,
+                                         const ModifierEvalContext *ctx,
+                                         Mesh *mesh)
+{
+  Mesh *result;
+  const SolidifyModifierData *smd = (SolidifyModifierData *)md;
+
+  MVert *mv, *mvert, *orig_mvert;
+  MEdge *ed, *medge, *orig_medge;
+  MLoop *ml, *mloop, *orig_mloop;
+  MPoly *mp, *mpoly, *orig_mpoly;
+  const uint numVerts = (uint)mesh->totvert;
+  const uint numEdges = (uint)mesh->totedge;
+  const uint numPolys = (uint)mesh->totpoly;
+  const uint numLoops = (uint)mesh->totloop;
+  uint newLoops = 0, newPolys = 0, newEdges = 0, newVerts = 0, rimVerts = 0;
+
+  /* only use material offsets if we have 2 or more materials  */
+  const short mat_nr_max = ctx->object->totcol > 1 ? ctx->object->totcol - 1 : 0;
+  const short mat_ofs = mat_nr_max ? smd->mat_ofs : 0;
+  const short mat_ofs_rim = mat_nr_max ? smd->mat_ofs_rim : 0;
+
+  /* use for edges */
+  /* over-alloc new_vert_arr, old_vert_arr */
+  uint *new_vert_arr = NULL;
+  STACK_DECLARE(new_vert_arr);
+
+  uint *new_edge_arr = NULL;
+  STACK_DECLARE(new_edge_arr);
+
+  uint *old_vert_arr = MEM_calloc_arrayN(
+      numVerts, sizeof(*old_vert_arr), "old_vert_arr in solidify");
+
+  uint *edge_users = NULL;
+  char *edge_order = NULL;
+
+  float(*vert_nors)[3] = NULL;
+  float(*poly_nors)[3] = NULL;
+
+  const bool need_poly_normals = (smd->flag & MOD_SOLIDIFY_NORMAL_CALC) ||
+                                 (smd->flag & MOD_SOLIDIFY_EVEN) ||
+                                 (smd->flag & MOD_SOLIDIFY_OFFSET_ANGLE_CLAMP);
+
+  const float ofs_orig = -(((-smd->offset_fac + 1.0f) * 0.5f) * smd->offset);
+  const float ofs_new = smd->offset + ofs_orig;
+  const float offset_fac_vg = smd->offset_fac_vg;
+  const float offset_fac_vg_inv = 1.0f - smd->offset_fac_vg;
+  const bool do_flip = (smd->flag & MOD_SOLIDIFY_FLIP) != 0;
+  const bool do_clamp = (smd->offset_clamp != 0.0f);
+  const bool do_angle_clamp = (smd->flag & MOD_SOLIDIFY_OFFSET_ANGLE_CLAMP) != 0;
+  const bool do_shell = ((smd->flag & MOD_SOLIDIFY_RIM) && (smd->flag & MOD_SOLIDIFY_NOSHELL)) ==
+                        0;
+
+  /* weights */
+  MDeformVert *dvert;
+  const bool defgrp_invert = (smd->flag & MOD_SOLIDIFY_VGROUP_INV) != 0;
+  int defgrp_index;
+
+  /* array size is doubled in case of using a shell */
+  const uint stride = do_shell ? 2 : 1;
+
+  MOD_get_vgroup(ctx->object, mesh, smd->defgrp_name, &dvert, &defgrp_index);
+
+  orig_mvert = mesh->mvert;
+  orig_medge = mesh->medge;
+  orig_mloop = mesh->mloop;
+  orig_mpoly = mesh->mpoly;
+
+  if (need_poly_normals) {
+    /* calculate only face normals */
+    poly_nors = MEM_malloc_arrayN(numPolys, sizeof(*poly_nors), __func__);
+    BKE_mesh_calc_normals_poly(orig_mvert,
+                               NULL,
+                               (int)numVerts,
+                               orig_mloop,
+                               orig_mpoly,
+                               (int)numLoops,
+                               (int)numPolys,
+                               poly_nors,
+                               true);
+  }
+
+  STACK_INIT(new_vert_arr, numVerts * 2);
+  STACK_INIT(new_edge_arr, numEdges * 2);
+
+  if (smd->flag & MOD_SOLIDIFY_RIM) {
+    BLI_bitmap *orig_mvert_tag = BLI_BITMAP_NEW(numVerts, __func__);
+    uint eidx;
+    uint i;
+
+#define INVALID_UNUSED ((uint)-1)
+#define INVALID_PAIR ((uint)-2)
+
+    new_vert_arr = MEM_malloc_arrayN(numVerts, 2 * sizeof(*new_vert_arr), __func__);
+    new_edge_arr = MEM_malloc_arrayN(((numEdges * 2) + numVerts), sizeof(*new_edge_arr), __func__);
+
+    edge_users = MEM_malloc_arrayN(numEdges, sizeof(*edge_users), "solid_mod edges");
+    edge_order = MEM_malloc_arrayN(numEdges, sizeof(*edge_order), "solid_mod eorder");
+
+    /* save doing 2 loops here... */
+#if 0
+    copy_vn_i(edge_users, numEdges, INVALID_UNUSED);
+#endif
+
+    for (eidx = 0, ed = orig_medge; eidx < numEdges; eidx++, ed++) {
+      edge_users[eidx] = INVALID_UNUSED;
+    }
+
+    for (i = 0, mp = orig_mpoly; i < numPolys; i++, mp++) {
+      MLoop *ml_prev;
+      int j;
+
+      ml = orig_mloop + mp->loopstart;
+      ml_prev = ml + (mp->totloop - 1);
+
+      for (j = 0; j < mp->totloop; j++, ml++) {
+        /* add edge user */
+        eidx = ml_prev->e;
+        if (edge_users[eidx] == INVALID_UNUSED) {
+          ed = orig_medge + eidx;
+          BLI_assert(ELEM(ml_prev->v, ed->v1, ed->v2) && ELEM(ml->v, ed->v1, ed->v2));
+          edge_users[eidx] = (ml_prev->v > ml->v) == (ed->v1 < ed->v2) ? i : (i + numPolys);
+          edge_order[eidx] = j;
+        }
+        else {
+          edge_users[eidx] = INVALID_PAIR;
+        }
+        ml_prev = ml;
+      }
+    }
+
+    for (eidx = 0, ed = orig_medge; eidx < numEdges; eidx++, ed++) {
+      if (!ELEM(edge_users[eidx], INVALID_UNUSED, INVALID_PAIR)) {
+        BLI_BITMAP_ENABLE(orig_mvert_tag, ed->v1);
+        BLI_BITMAP_ENABLE(orig_mvert_tag, ed->v2);
+        STACK_PUSH(new_edge_arr, eidx);
+        newPolys++;
+        newLoops += 4;
+      }
+    }
+
+    for (i = 0; i < numVerts; i++) {
+      if (BLI_BITMAP_TEST(orig_mvert_tag, i)) {
+        old_vert_arr[i] = STACK_SIZE(new_vert_arr);
+        STACK_PUSH(new_vert_arr, i);
+        rimVerts++;
+      }
+      else {
+        old_vert_arr[i] = INVALID_UNUSED;
+      }
+    }
+
+    MEM_freeN(orig_mvert_tag);
+  }
+
+  if (do_shell == false) {
+    /* only add rim vertices */
+    newVerts = rimVerts;
+    /* each extruded face needs an opposite edge */
+    newEdges = newPolys;
+  }
+  else {
+    /* (stride == 2) in this case, so no need to add newVerts/newEdges */
+    BLI_assert(newVerts == 0);
+    BLI_assert(newEdges == 0);
+  }
+
+  if (smd->flag & MOD_SOLIDIFY_NORMAL_CALC) {
+    vert_nors = MEM_calloc_arrayN(numVerts, 3 * sizeof(float), "mod_solid_vno_hq");
+    mesh_calc_hq_normal(mesh, poly_nors, vert_nors);
+  }
+
+  result = BKE_mesh_new_nomain_from_template(mesh,
+                                             (int)((numVerts * stride) + newVerts),
+                                             (int)((numEdges * stride) + newEdges + rimVerts),
+                                             0,
+                                             (int)((numLoops * stride) + newLoops),
+                                             (int)((numPolys * stride) + newPolys));
+
+  mpoly = result->mpoly;
+  mloop = result->mloop;
+  medge = result->medge;
+  mvert = result->mvert;
+
+  if (do_shell) {
+    CustomData_copy_data(&mesh->vdata, &result->vdata, 0, 0, (int)numVerts);
+    CustomData_copy_data(&mesh->vdata, &result->vdata, 0, (int)numVerts, (int)numVerts);
+
+    CustomData_copy_data(&mesh->edata, &result->edata, 0, 0, (int)numEdges);
+    CustomData_copy_data(&mesh->edata, &result->edata, 0, (int)numEdges, (int)numEdges);
+
+    CustomData_copy_data(&mesh->ldata, &result->ldata, 0, 0, (int)numLoops);
+    /* DO NOT copy here the 'copied' part of loop data, we want to reverse loops
+     * (so that winding of copied face get reversed, so that normals get reversed
+     * and point in expected direction...).
+     * If we also copy data here, then this data get overwritten
+     * (and allocated memory becomes memleak). */
+
+    CustomData_copy_data(&mesh->pdata, &result->pdata, 0, 0, (int)numPolys);
+    CustomData_copy_data(&mesh->pdata, &result->pdata, 0, (int)numPolys, (int)numPolys);
+  }
+  else {
+    int i, j;
+    CustomData_copy_data(&mesh->vdata, &result->vdata, 0, 0, (int)numVerts);
+    for (i = 0, j = (int)numVerts; i < numVerts; i++) {
+      if (old_vert_arr[i] != INVALID_UNUSED) {
+        CustomData_copy_data(&mesh->vdata, &result->vdata, i, j, 1);
+        j++;
+      }
+    }
+
+    CustomData_copy_data(&mesh->edata, &result->edata, 0, 0, (int)numEdges);
+
+    for (i = 0, j = (int)numEdges; i < numEdges; i++) {
+      if (!ELEM(edge_users[i], INVALID_UNUSED, INVALID_PAIR)) {
+        MEdge *ed_src, *ed_dst;
+        CustomData_copy_data(&mesh->edata, &result->edata, i, j, 1);
+
+        ed_src = &medge[i];
+        ed_dst = &medge[j];
+        ed_dst->v1 = old_vert_arr[ed_src->v1] + numVerts;
+        ed_dst->v2 = old_vert_arr[ed_src->v2] + numVerts;
+        j++;
+      }
+    }
+
+    /* will be created later */
+    CustomData_copy_data(&mesh->ldata, &result->ldata, 0, 0, (int)numLoops);
+    CustomData_copy_data(&mesh->pdata, &result->pdata, 0, 0, (int)numPolys);
+  }
+
+  /* initializes: (i_end, do_shell_align, mv)  */
+#define INIT_VERT_ARRAY_OFFSETS(test) \
+  if (((ofs_new >= ofs_orig) == do_flip) == test) { \
+    i_end = numVerts; \
+    do_shell_align = true; \
+    mv = mvert; \
+  } \
+  else { \
+    if (do_shell) { \
+      i_end = numVerts; \
+      do_shell_align = true; \
+    } \
+    else { \
+      i_end = newVerts; \
+      do_shell_align = false; \
+    } \
+    mv = &mvert[numVerts]; \
+  } \
+  (void)0
+
+  /* flip normals */
+
+  if (do_shell) {
+    uint i;
+
+    mp = mpoly + numPolys;
+    for (i = 0; i < mesh->totpoly; i++, mp++) {
+      const int loop_end = mp->totloop - 1;
+      MLoop *ml2;
+      uint e;
+      int j;
+
+      /* reverses the loop direction (MLoop.v as well as custom-data)
+       * MLoop.e also needs to be corrected too, done in a separate loop below. */
+      ml2 = mloop + mp->loopstart + mesh->totloop;
+#if 0
+      for (j = 0; j < mp->totloop; j++) {
+        CustomData_copy_data(&mesh->ldata,
+                             &result->ldata,
+                             mp->loopstart + j,
+                             mp->loopstart + (loop_end - j) + mesh->totloop,
+                             1);
+      }
+#else
+      /* slightly more involved, keep the first vertex the same for the copy,
+       * ensures the diagonals in the new face match the original. */
+      j = 0;
+      for (int j_prev = loop_end; j < mp->totloop; j_prev = j++) {
+        CustomData_copy_data(&mesh->ldata,
+                             &result->ldata,
+                             mp->loopstart + j,
+                             mp->loopstart + (loop_end - j_prev) + mesh->totloop,
+                             1);
+      }
+#endif
+
+      if (mat_ofs) {
+        mp->mat_nr += mat_ofs;
+        CLAMP(mp->mat_nr, 0, mat_nr_max);
+      }
+
+      e = ml2[0].e;
+      for (j = 0; j < loop_end; j++) {
+        ml2[j].e = ml2[j + 1].e;
+      }
+      ml2[loop_end].e = e;
+
+      mp->loopstart += mesh->totloop;
+
+      for (j = 0; j < mp->totloop; j++) {
+        ml2[j].e += numEdges;
+        ml2[j].v += numVerts;
+      }
+    }
+
+    for (i = 0, ed = medge + numEdges; i < numEdges; i++, ed++) {
+      ed->v1 += numVerts;
+      ed->v2 += numVerts;
+    }
+  }
+
+  /* note, copied vertex layers don't have flipped normals yet. do this after applying offset */
+  if ((smd->flag & MOD_SOLIDIFY_EVEN) == 0) {
+    /* no even thickness, very simple */
+    float scalar_short;
+    float scalar_short_vgroup;
+
+    /* for clamping */
+    float *vert_lens = NULL;
+    float *vert_angs = NULL;
+    const float offset = fabsf(smd->offset) * smd->offset_clamp;
+    const float offset_sq = offset * offset;
+
+    if (do_clamp) {
+      uint i;
+
+      vert_lens = MEM_malloc_arrayN(numVerts, sizeof(float), "vert_lens");
+      copy_vn_fl(vert_lens, (int)numVerts, FLT_MAX);
+      for (i = 0; i < numEdges; i++) {
+        const float ed_len_sq = len_squared_v3v3(mvert[medge[i].v1].co, mvert[medge[i].v2].co);
+        vert_lens[medge[i].v1] = min_ff(vert_lens[medge[i].v1], ed_len_sq);
+        vert_lens[medge[i].v2] = min_ff(vert_lens[medge[i].v2], ed_len_sq);
+      }
+      if (do_angle_clamp) {
+        uint eidx;
+        vert_angs = MEM_malloc_arrayN(numVerts, sizeof(float), "vert_angs");
+        copy_vn_fl(vert_angs, (int)numVerts, 0.5f * M_PI);
+        uint(*edge_user_pairs)[2] = MEM_malloc_arrayN(
+            numEdges, sizeof(*edge_user_pairs), "edge_user_pairs");
+        for (eidx = 0; eidx < numEdges; eidx++) {
+          edge_user_pairs[eidx][0] = INVALID_UNUSED;
+          edge_user_pairs[eidx][1] = INVALID_UNUSED;
+        }
+        for (i = 0, mp = orig_mpoly; i < numPolys; i++, mp++) {
+          ml = orig_mloop + mp->loopstart;
+          MLoop *ml_prev = ml + (mp->totloop - 1);
+
+          for (int j = 0; j < mp->totloop; j++, ml++) {
+            /* add edge user */
+            eidx = ml_prev->e;
+            ed = orig_medge + eidx;
+            BLI_assert(ELEM(ml_prev->v, ed->v1, ed->v2) && ELEM(ml->v, ed->v1, ed->v2));
+            char flip = (char)((ml_prev->v > ml->v) == (ed->v1 < ed->v2));
+            if (edge_user_pairs[eidx][flip] == INVALID_UNUSED) {
+              edge_user_pairs[eidx][flip] = i;
+            }
+            else {
+              edge_user_pairs[eidx][0] = INVALID_PAIR;
+              edge_user_pairs[eidx][1] = INVALID_PAIR;
+            }
+            ml_prev = ml;
+          }
+        }
+        ed = orig_medge;
+        float e[3];
+        for (i = 0; i < numEdges; i++, ed++) {
+          if (!ELEM(edge_user_pairs[i][0], INVALID_UNUSED, INVALID_PAIR) &&
+              !ELEM(edge_user_pairs[i][1], INVALID_UNUSED, INVALID_PAIR)) {
+            const float *n0 = poly_nors[edge_user_pairs[i][0]];
+            const float *n1 = poly_nors[edge_user_pairs[i][1]];
+            sub_v3_v3v3(e, orig_mvert[ed->v1].co, orig_mvert[ed->v2].co);
+            normalize_v3(e);
+            const float angle = angle_signed_on_axis_v3v3_v3(n0, n1, e);
+            vert_angs[ed->v1] = max_ff(vert_angs[ed->v1], angle);
+            vert_angs[ed->v2] = max_ff(vert_angs[ed->v2], angle);
+          }
+        }
+        MEM_freeN(edge_user_pairs);
+      }
+    }
+
+    if (ofs_new != 0.0f) {
+      uint i_orig, i_end;
+      bool do_shell_align;
+
+      scalar_short = scalar_short_vgroup = ofs_new / 32767.0f;
+
+      INIT_VERT_ARRAY_OFFSETS(false);
+
+      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
+        const uint i = do_shell_align ? i_orig : new_vert_arr[i_orig];
+        if (dvert) {
+          MDeformVert *dv = &dvert[i];
+          if (defgrp_invert) {
+            scalar_short_vgroup = 1.0f - defvert_find_weight(dv, defgrp_index);
+          }
+          else {
+            scalar_short_vgroup = defvert_find_weight(dv, defgrp_index);
+          }
+          scalar_short_vgroup = (offset_fac_vg + (scalar_short_vgroup * offset_fac_vg_inv)) *
+                                scalar_short;
+        }
+        if (do_clamp && offset > FLT_EPSILON) {
+          /* always reset because we may have set before */
+          if (dvert == NULL) {
+            scalar_short_vgroup = scalar_short;
+          }
+          if (do_angle_clamp) {
+            float cos_ang = cosf(((2 * M_PI) - vert_angs[i]) * 0.5f);
+            if (cos_ang > 0) {
+              float max_off = sqrtf(vert_lens[i]) * 0.5f / cos_ang;
+              if (max_off < offset * 0.5f) {
+                scalar_short_vgroup *= max_off / offset * 2;
+              }
+            }
+          }
+          else {
+            if (vert_lens[i] < offset_sq) {
+              float scalar = sqrtf(vert_lens[i]) / offset;
+              scalar_short_vgroup *= scalar;
+            }
+          }
+        }
+        madd_v3v3short_fl(mv->co, mv->no, scalar_short_vgroup);
+      }
+    }
+
+    if (ofs_orig != 0.0f) {
+      uint i_orig, i_end;
+      bool do_shell_align;
+
+      scalar_short = scalar_short_vgroup = ofs_orig / 32767.0f;
+
+      /* as above but swapped */
+      INIT_VERT_ARRAY_OFFSETS(true);
+
+      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
+        const uint i = do_shell_align ? i_orig : new_vert_arr[i_orig];
+        if (dvert) {
+          MDeformVert *dv = &dvert[i];
+          if (defgrp_invert) {
+            scalar_short_vgroup = 1.0f - defvert_find_weight(dv, defgrp_index);
+          }
+          else {
+            scalar_short_vgroup = defvert_find_weight(dv, defgrp_index);
+          }
+          scalar_short_vgroup = (offset_fac_vg + (scalar_short_vgroup * offset_fac_vg_inv)) *
+                                scalar_short;
+        }
+        if (do_clamp && offset > FLT_EPSILON) {
+          /* always reset because we may have set before */
+          if (dvert == NULL) {
+            scalar_short_vgroup = scalar_short;
+          }
+          if (do_angle_clamp) {
+            float cos_ang = cosf(vert_angs[i_orig] * 0.5f);
+            if (cos_ang > 0) {
+              float max_off = sqrtf(vert_lens[i]) * 0.5f / cos_ang;
+              if (max_off < offset * 0.5f) {
+                scalar_short_vgroup *= max_off / offset * 2;
+              }
+            }
+          }
+          else {
+            if (vert_lens[i] < offset_sq) {
+              float scalar = sqrtf(vert_lens[i]) / offset;
+              scalar_short_vgroup *= scalar;
+            }
+          }
+        }
+        madd_v3v3short_fl(mv->co, mv->no, scalar_short_vgroup);
+      }
+    }
+
+    if (do_clamp) {
+      MEM_freeN(vert_lens);
+      if (do_angle_clamp) {
+        MEM_freeN(vert_angs);
+      }
+    }
+  }
+  else {
+#ifdef USE_NONMANIFOLD_WORKAROUND
+    const bool check_non_manifold = (smd->flag & MOD_SOLIDIFY_NORMAL_CALC) != 0;
+#endif
+    /* same as EM_solidify() in editmesh_lib.c */
+    float *vert_angles = MEM_calloc_arrayN(
+        numVerts, 2 * sizeof(float), "mod_solid_pair"); /* 2 in 1 */
+    float *vert_accum = vert_angles + numVerts;
+    uint vidx;
+    uint i;
+
+    if (vert_nors == NULL) {
+      vert_nors = MEM_malloc_arrayN(numVerts, 3 * sizeof(float), "mod_solid_vno");
+      for (i = 0, mv = mvert; i < numVerts; i++, mv++) {
+        normal_short_to_float_v3(vert_nors[i], mv->no);
+      }
+    }
+
+    for (i = 0, mp = mpoly; i < numPolys; i++, mp++) {
+      /* #BKE_mesh_calc_poly_angles logic is inlined here */
+      float nor_prev[3];
+      float nor_next[3];
+
+      int i_curr = mp->totloop - 1;
+      int i_next = 0;
+
+      ml = &mloop[mp->loopstart];
+
+      sub_v3_v3v3(nor_prev, mvert[ml[i_curr - 1].v].co, mvert[ml[i_curr].v].co);
+      normalize_v3(nor_prev);
+
+      while (i_next < mp->totloop) {
+        float angle;
+        sub_v3_v3v3(nor_next, mvert[ml[i_curr].v].co, mvert[ml[i_next].v].co);
+        normalize_v3(nor_next);
+        angle = angle_normalized_v3v3(nor_prev, nor_next);
+
+        /* --- not related to angle calc --- */
+        if (angle < FLT_EPSILON) {
+          angle = FLT_EPSILON;
+        }
+
+        vidx = ml[i_curr].v;
+        vert_accum[vidx] += angle;
+
+#ifdef USE_NONMANIFOLD_WORKAROUND
+        /* skip 3+ face user edges */
+        if ((check_non_manifold == false) ||
+            LIKELY(((orig_medge[ml[i_curr].e].flag & ME_EDGE_TMP_TAG) == 0) &&
+                   ((orig_medge[ml[i_next].e].flag & ME_EDGE_TMP_TAG) == 0))) {
+          vert_angles[vidx] += shell_v3v3_normalized_to_dist(vert_nors[vidx], poly_nors[i]) *
+                               angle;
+        }
+        else {
+          vert_angles[vidx] += angle;
+        }
+#else
+        vert_angles[vidx] += shell_v3v3_normalized_to_dist(vert_nors[vidx], poly_nors[i]) * angle;
+#endif
+        /* --- end non-angle-calc section --- */
+
+        /* step */
+        copy_v3_v3(nor_prev, nor_next);
+        i_curr = i_next;
+        i_next++;
+      }
+    }
+
+    /* vertex group support */
+    if (dvert) {
+      MDeformVert *dv = dvert;
+      float scalar;
+
+      if (defgrp_invert) {
+        for (i = 0; i < numVerts; i++, dv++) {
+          scalar = 1.0f - defvert_find_weight(dv, defgrp_index);
+          scalar = offset_fac_vg + (scalar * offset_fac_vg_inv);
+          vert_angles[i] *= scalar;
+        }
+      }
+      else {
+        for (i = 0; i < numVerts; i++, dv++) {
+          scalar = defvert_find_weight(dv, defgrp_index);
+          scalar = offset_fac_vg + (scalar * offset_fac_vg_inv);
+          vert_angles[i] *= scalar;
+        }
+      }
+    }
+
+    if (do_clamp) {
+      const float clamp_fac = 1 + (do_angle_clamp ? fabsf(smd->offset_fac) : 0);
+      const float offset = fabsf(smd->offset) * smd->offset_clamp * clamp_fac;
+      if (offset > FLT_EPSILON) {
+        float *vert_lens_sq = MEM_malloc_arrayN(numVerts, sizeof(float), "vert_lens_sq");
+        const float offset_sq = offset * offset;
+        copy_vn_fl(vert_lens_sq, (int)numVerts, FLT_MAX);
+        for (i = 0; i < numEdges; i++) {
+          const float ed_len = len_squared_v3v3(mvert[medge[i].v1].co, mvert[medge[i].v2].co);
+          vert_lens_sq[medge[i].v1] = min_ff(vert_lens_sq[medge[i].v1], ed_len);
+          vert_lens_sq[medge[i].v2] = min_ff(vert_lens_sq[medge[i].v2], ed_len);
+        }
+        if (do_angle_clamp) {
+          uint eidx;
+          float *vert_angs = MEM_malloc_arrayN(numVerts, sizeof(float), "vert_angs even");
+          copy_vn_fl(vert_angs, (int)numVerts, 0.5f * M_PI);
+          uint(*edge_user_pairs)[2] = MEM_malloc_arrayN(
+              numEdges, sizeof(*edge_user_pairs), "edge_user_pairs");
+          for (eidx = 0; eidx < numEdges; eidx++) {
+            edge_user_pairs[eidx][0] = INVALID_UNUSED;
+            edge_user_pairs[eidx][1] = INVALID_UNUSED;
+          }
+          for (i = 0, mp = orig_mpoly; i < numPolys; i++, mp++) {
+            ml = orig_mloop + mp->loopstart;
+            MLoop *ml_prev = ml + (mp->totloop - 1);
+
+            for (int j = 0; j < mp->totloop; j++, ml++) {
+              /* add edge user */
+              eidx = ml_prev->e;
+              ed = orig_medge + eidx;
+              BLI_assert(ELEM(ml_prev->v, ed->v1, ed->v2) && ELEM(ml->v, ed->v1, ed->v2));
+              char flip = (char)((ml_prev->v > ml->v) == (ed->v1 < ed->v2));
+              if (edge_user_pairs[eidx][flip] == INVALID_UNUSED) {
+                edge_user_pairs[eidx][flip] = i;
+              }
+              else {
+                edge_user_pairs[eidx][0] = INVALID_PAIR;
+                edge_user_pairs[eidx][1] = INVALID_PAIR;
+              }
+              ml_prev = ml;
+            }
+          }
+          ed = orig_medge;
+          for (i = 0; i < numEdges; i++, ed++) {
+            if (!ELEM(edge_user_pairs[i][0], INVALID_UNUSED, INVALID_PAIR) &&
+                !ELEM(edge_user_pairs[i][1], INVALID_UNUSED, INVALID_PAIR)) {
+              const float *n0 = poly_nors[edge_user_pairs[i][0]];
+              const float *n1 = poly_nors[edge_user_pairs[i][1]];
+              const float angle = M_PI - angle_normalized_v3v3(n0, n1);
+              vert_angs[ed->v1] = max_ff(vert_angs[ed->v1], angle);
+              vert_angs[ed->v2] = max_ff(vert_angs[ed->v2], angle);
+            }
+          }
+          MEM_freeN(edge_user_pairs);
+
+          for (i = 0; i < numVerts; i++) {
+            float cos_ang = cosf(vert_angs[i] * 0.5f);
+            if (cos_ang > 0) {
+              float max_off = sqrtf(vert_lens_sq[i]) * 0.5f / cos_ang;
+              if (max_off < offset * 0.5f) {
+                vert_angles[i] *= max_off / offset * 2;
+              }
+            }
+          }
+          MEM_freeN(vert_angs);
+        }
+        else {
+          for (i = 0; i < numVerts; i++) {
+            if (vert_lens_sq[i] < offset_sq) {
+              float scalar = sqrtf(vert_lens_sq[i]) / offset;
+              vert_angles[i] *= scalar;
+            }
+          }
+        }
+        MEM_freeN(vert_lens_sq);
+      }
+    }
+
+#undef INVALID_UNUSED
+#undef INVALID_PAIR
+
+    if (ofs_new != 0.0f) {
+      uint i_orig, i_end;
+      bool do_shell_align;
+
+      INIT_VERT_ARRAY_OFFSETS(false);
+
+      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
+        const uint i_other = do_shell_align ? i_orig : new_vert_arr[i_orig];
+        if (vert_accum[i_other]) { /* zero if unselected */
+          madd_v3_v3fl(
+              mv->co, vert_nors[i_other], ofs_new * (vert_angles[i_other] / vert_accum[i_other]));
+        }
+      }
+    }
+
+    if (ofs_orig != 0.0f) {
+      uint i_orig, i_end;
+      bool do_shell_align;
+
+      /* same as above but swapped, intentional use of 'ofs_new' */
+      INIT_VERT_ARRAY_OFFSETS(true);
+
+      for (i_orig = 0; i_orig < i_end; i_orig++, mv++) {
+        const uint i_other = do_shell_align ? i_orig : new_vert_arr[i_orig];
+        if (vert_accum[i_other]) { /* zero if unselected */
+          madd_v3_v3fl(
+              mv->co, vert_nors[i_other], ofs_orig * (vert_angles[i_other] / vert_accum[i_other]));
+        }
+      }
+    }
+
+    MEM_freeN(vert_angles);
+  }
+
+  if (vert_nors) {
+    MEM_freeN(vert_nors);
+  }
+
+  /* must recalculate normals with vgroups since they can displace unevenly [#26888] */
+  if ((mesh->runtime.cd_dirty_vert & CD_MASK_NORMAL) || (smd->flag & MOD_SOLIDIFY_RIM) || dvert) {
+    result->runtime.cd_dirty_vert |= CD_MASK_NORMAL;
+  }
+  else if (do_shell) {
+    uint i;
+    /* flip vertex normals for copied verts */
+    mv = mvert + numVerts;
+    for (i = 0; i < numVerts; i++, mv++) {
+      negate_v3_short(mv->no);
+    }
+  }
+
+  if (smd->flag & MOD_SOLIDIFY_RIM) {
+    uint i;
+
+    /* bugger, need to re-calculate the normals for the new edge faces.
+     * This could be done in many ways, but probably the quickest way
+     * is to calculate the average normals for side faces only.
+     * Then blend them with the normals of the edge verts.
+     *
+     * at the moment its easiest to allocate an entire array for every vertex,
+     * even though we only need edge verts - campbell
+     */
+
+#define SOLIDIFY_SIDE_NORMALS
+
+#ifdef SOLIDIFY_SIDE_NORMALS
+    /* Note that, due to the code setting cd_dirty_vert a few lines above,
+     * do_side_normals is always false. - Sybren */
+    const bool do_side_normals = !(result->runtime.cd_dirty_vert & CD_MASK_NORMAL);
+    /* annoying to allocate these since we only need the edge verts, */
+    float(*edge_vert_nos)[3] = do_side_normals ?
+                                   MEM_calloc_arrayN(numVerts, 3 * sizeof(float), __func__) :
+                                   NULL;
+    float nor[3];
+#endif
+    const uchar crease_rim = smd->crease_rim * 255.0f;
+    const uchar crease_outer = smd->crease_outer * 255.0f;
+    const uchar crease_inner = smd->crease_inner * 255.0f;
+
+    int *origindex_edge;
+    int *orig_ed;
+    uint j;
+
+    if (crease_rim || crease_outer || crease_inner) {
+      result->cd_flag |= ME_CDFLAG_EDGE_CREASE;
+    }
+
+    /* add faces & edges */
+    origindex_edge = CustomData_get_layer(&result->edata, CD_ORIGINDEX);
+    orig_ed = (origindex_edge) ? &origindex_edge[(numEdges * stride) + newEdges] : NULL;
+    ed = &medge[(numEdges * stride) + newEdges]; /* start after copied edges */
+    for (i = 0; i < rimVerts; i++, ed++) {
+      ed->v1 = new_vert_arr[i];
+      ed->v2 = (do_shell ? new_vert_arr[i] : i) + numVerts;
+      ed->flag |= ME_EDGEDRAW | ME_EDGERENDER;
+
+      if (orig_ed) {
+        *orig_ed = ORIGINDEX_NONE;
+        orig_ed++;
+      }
+
+      if (crease_rim) {
+        ed->crease = crease_rim;
+      }
+    }
+
+    /* faces */
+    mp = mpoly + (numPolys * stride);
+    ml = mloop + (numLoops * stride);
+    j = 0;
+    for (i = 0; i < newPolys; i++, mp++) {
+      uint eidx = new_edge_arr[i];
+      uint pidx = edge_users[eidx];
+      int k1, k2;
+      bool flip;
+
+      if (pidx >= numPolys) {
+        pidx -= numPolys;
+        flip = true;
+      }
+      else {
+        flip = false;
+      }
+
+      ed = medge + eidx;
+
+      /* copy most of the face settings */
+      CustomData_copy_data(
+          &mesh->pdata, &result->pdata, (int)pidx, (int)((numPolys * stride) + i), 1);
+      mp->loopstart = (int)(j + (numLoops * stride));
+      mp->flag = mpoly[pidx].flag;
+
+      /* notice we use 'mp->totloop' which is later overwritten,
+       * we could lookup the original face but there's no point since this is a copy
+       * and will have the same value, just take care when changing order of assignment */
+
+      /* prev loop */
+      k1 = mpoly[pidx].loopstart + (((edge_order[eidx] - 1) + mp->totloop) % mp->totloop);
+
+      k2 = mpoly[pidx].loopstart + (edge_order[eidx]);
+
+      mp->totloop = 4;
+
+      CustomData_copy_data(
+          &mesh->ldata, &result->ldata, k2, (int)((numLoops * stride) + j + 0), 1);
+      CustomData_copy_data(
+          &mesh->ldata, &result->ldata, k1, (int)((numLoops * stride) + j + 1), 1);
+      CustomData_copy_data(
+          &mesh->ldata, &result->ldata, k1, (int)((numLoops * stride) + j + 2), 1);
+      CustomData_copy_data(
+          &mesh->ldata, &result->ldata, k2, (int)((numLoops * stride) + j + 3), 1);
+
+      if (flip == false) {
+        ml[j].v = ed->v1;
+        ml[j++].e = eidx;
+
+        ml[j].v = ed->v2;
+        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v2] + newEdges;
+
+        ml[j].v = (do_shell ? ed->v2 : old_vert_arr[ed->v2]) + numVerts;
+        ml[j++].e = (do_shell ? eidx : i) + numEdges;
+
+        ml[j].v = (do_shell ? ed->v1 : old_vert_arr[ed->v1]) + numVerts;
+        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v1] + newEdges;
+      }
+      else {
+        ml[j].v = ed->v2;
+        ml[j++].e = eidx;
+
+        ml[j].v = ed->v1;
+        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v1] + newEdges;
+
+        ml[j].v = (do_shell ? ed->v1 : old_vert_arr[ed->v1]) + numVerts;
+        ml[j++].e = (do_shell ? eidx : i) + numEdges;
+
+        ml[j].v = (do_shell ? ed->v2 : old_vert_arr[ed->v2]) + numVerts;
+        ml[j++].e = (numEdges * stride) + old_vert_arr[ed->v2] + newEdges;
+      }
+
+      if (origindex_edge) {
+        origindex_edge[ml[j - 3].e] = ORIGINDEX_NONE;
+        origindex_edge[ml[j - 1].e] = ORIGINDEX_NONE;
+      }
+
+      /* use the next material index if option enabled */
+      if (mat_ofs_rim) {
+        mp->mat_nr += mat_ofs_rim;
+        CLAMP(mp->mat_nr, 0, mat_nr_max);
+      }
+      if (crease_outer) {
+        /* crease += crease_outer; without wrapping */
+        char *cr = &(ed->crease);
+        int tcr = *cr + crease_outer;
+        *cr = tcr > 255 ? 255 : tcr;
+      }
+
+      if (crease_inner) {
+        /* crease += crease_inner; without wrapping */
+        char *cr = &(medge[numEdges + (do_shell ? eidx : i)].crease);
+        int tcr = *cr + crease_inner;
+        *cr = tcr > 255 ? 255 : tcr;
+      }
+
+#ifdef SOLIDIFY_SIDE_NORMALS
+      if (do_side_normals) {
+        normal_quad_v3(nor,
+                       mvert[ml[j - 4].v].co,
+                       mvert[ml[j - 3].v].co,
+                       mvert[ml[j - 2].v].co,
+                       mvert[ml[j - 1].v].co);
+
+        add_v3_v3(edge_vert_nos[ed->v1], nor);
+        add_v3_v3(edge_vert_nos[ed->v2], nor);
+      }
+#endif
+    }
+
+#ifdef SOLIDIFY_SIDE_NORMALS
+    if (do_side_normals) {
+      const MEdge *ed_orig = medge;
+      ed = medge + (numEdges * stride);
+      for (i = 0; i < rimVerts; i++, ed++, ed_orig++) {
+        float nor_cpy[3];
+        short *nor_short;
+        int k;
+
+        /* note, only the first vertex (lower half of the index) is calculated */
+        BLI_assert(ed->v1 < numVerts);
+        normalize_v3_v3(nor_cpy, edge_vert_nos[ed_orig->v1]);
+
+        for (k = 0; k < 2; k++) { /* loop over both verts of the edge */
+          nor_short = mvert[*(&ed->v1 + k)].no;
+          normal_short_to_float_v3(nor, nor_short);
+          add_v3_v3(nor, nor_cpy);
+          normalize_v3(nor);
+          normal_float_to_short_v3(nor_short, nor);
+        }
+      }
+
+      MEM_freeN(edge_vert_nos);
+    }
+#endif
+
+    MEM_freeN(new_vert_arr);
+    MEM_freeN(new_edge_arr);
+
+    MEM_freeN(edge_users);
+    MEM_freeN(edge_order);
+  }
+
+  if (old_vert_arr) {
+    MEM_freeN(old_vert_arr);
+  }
+
+  if (poly_nors) {
+    MEM_freeN(poly_nors);
+  }
+
+  if (numPolys == 0 && numVerts != 0) {
+    modifier_setError(md, "Faces needed for useful output");
+  }
+
+  return result;
+}
+
+#undef SOLIDIFY_SIDE_NORMALS
+
+/** \} */
diff --git a/source/blender/modifiers/intern/MOD_solidify_nonmanifold.c b/source/blender/modifiers/intern/MOD_solidify_nonmanifold.c
new file mode 100644
index 00000000000..be7bbb86e4d
--- /dev/null
+++ b/source/blender/modifiers/intern/MOD_solidify_nonmanifold.c
@@ -0,0 +1,2268 @@
+/*
+ * 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.
+ */
+
+/** \file
+ * \ingroup modifiers
+ */
+
+#include "BLI_utildefines.h"
+
+#include "BLI_math.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "BKE_mesh.h"
+#include "BKE_particle.h"
+#include "BKE_deform.h"
+
+#include "MOD_modifiertypes.h"
+#include "MOD_util.h"
+#include "MOD_solidify_util.h" /* Own include. */
+
+#ifdef __GNUC__
+#  pragma GCC diagnostic error "-Wsign-conversion"
+#endif
+
+/* -------------------------------------------------------------------- */
+/** \name Local Utilities
+ * \{ */
+
+/**
+ * Similar to #project_v3_v3v3_normalized that returns the dot-product.
+ */
+static float project_v3_v3(float r[3], const float a[3])
+{
+  float d = dot_v3v3(r, a);
+  r[0] -= a[0] * d;
+  r[1] -= a[1] * d;
+  r[2] -= a[2] * d;
+  return d;
+}
+
+static float angle_signed_on_axis_normalized_v3v3_v3(const float n[3],
+                                                     const float ref_n[3],
+                                                     const float axis[3])
+{
+  float d = dot_v3v3(n, ref_n);
+  CLAMP(d, -1, 1);
+  float angle = acosf(d);
+  float cross[3];
+  cross_v3_v3v3(cross, n, ref_n);
+  if (dot_v3v3(cross, axis) >= 0) {
+    angle = 2 * M_PI - angle;
+  }
+  return angle;
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Main Solidify Function
+ * \{ */
+
+/* Data structures for manifold solidify. */
+
+typedef struct NewFaceRef {
+  MPoly *face;
+  uint index;
+  bool reversed;
+  struct NewEdgeRef **link_edges;
+} NewFaceRef;
+
+typedef struct OldEdgeFaceRef {
+  uint *faces;
+  uint faces_len;
+  bool *faces_reversed;
+  uint used;
+} OldEdgeFaceRef;
+
+typedef struct OldVertEdgeRef {
+  uint *edges;
+  uint edges_len;
+} OldVertEdgeRef;
+
+typedef struct NewEdgeRef {
+  uint old_edge;
+  NewFaceRef *faces[2];
+  struct EdgeGroup *link_edge_groups[2];
+  float angle;
+  uint new_edge;
+} NewEdgeRef;
+
+typedef struct EdgeGroup {
+  bool valid;
+  NewEdgeRef **edges;
+  uint edges_len;
+  uint open_face_edge;
+  bool is_orig_closed;
+  bool is_even_split;
+  uint split;
+  bool is_singularity;
+  uint topo_group;
+  float co[3];
+  float no[3];
+  uint new_vert;
+} EdgeGroup;
+
+typedef struct FaceKeyPair {
+  float angle;
+  NewFaceRef *face;
+} FaceKeyPair;
+
+static int comp_float_int_pair(const void *a, const void *b)
+{
+  FaceKeyPair *x = (FaceKeyPair *)a;
+  FaceKeyPair *y = (FaceKeyPair *)b;
+  return (int)(x->angle > y->angle) - (int)(x->angle < y->angle);
+}
+
+Mesh *MOD_solidify_nonmanifold_applyModifier(ModifierData *md,
+                                             const ModifierEvalContext *ctx,
+                                             Mesh *mesh)
+{
+  Mesh *result;
+  const SolidifyModifierData *smd = (SolidifyModifierData *)md;
+
+  MVert *mv, *mvert, *orig_mvert;
+  MEdge *ed, *medge, *orig_medge;
+  MLoop *ml, *mloop, *orig_mloop;
+  MPoly *mp, *mpoly, *orig_mpoly;
+  const uint numVerts = (uint)mesh->totvert;
+  const uint numEdges = (uint)mesh->totedge;
+  const uint numPolys = (uint)mesh->totpoly;
+  const uint numLoops = (uint)mesh->totloop;
+
+  if (numPolys == 0 && numVerts != 0) {
+    modifier_setError(md, "Faces needed for useful output");
+    return mesh;
+  }
+
+  /* Only use material offsets if we have 2 or more materials. */
+  const short mat_nrs = ctx->object->totcol > 1 ? ctx->object->totcol : 1;
+  const short mat_nr_max = mat_nrs - 1;
+  const short mat_ofs = mat_nrs > 1 ? smd->mat_ofs : 0;
+  const short mat_ofs_rim = mat_nrs > 1 ? smd->mat_ofs_rim : 0;
+
+  float(*poly_nors)[3] = NULL;
+
+  const float ofs_front = (smd->offset_fac + 1.0f) * 0.5f * smd->offset;
+  const float ofs_back = ofs_front - smd->offset * smd->offset_fac;
+  const float offset_fac_vg = smd->offset_fac_vg;
+  const float offset_fac_vg_inv = 1.0f - smd->offset_fac_vg;
+  const float offset = fabsf(smd->offset) * smd->offset_clamp;
+  const bool do_angle_clamp = smd->flag & MOD_SOLIDIFY_OFFSET_ANGLE_CLAMP;
+  const bool do_flip = (smd->flag & MOD_SOLIDIFY_FLIP) != 0;
+  const bool do_rim = smd->flag & MOD_SOLIDIFY_RIM;
+  const bool do_shell = ((smd->flag & MOD_SOLIDIFY_RIM) && (smd->flag & MOD_SOLIDIFY_NOSHELL)) ==
+                        0;
+  const bool do_clamp = (smd->offset_clamp != 0.0f);
+
+  MDeformVert *dvert;
+  const bool defgrp_invert = (smd->flag & MOD_SOLIDIFY_VGROUP_INV) != 0;
+  int defgrp_index;
+
+  MOD_get_vgroup(ctx->object, mesh, smd->defgrp_name, &dvert, &defgrp_index);
+
+  orig_mvert = mesh->mvert;
+  orig_medge = mesh->medge;
+  orig_mloop = mesh->mloop;
+  orig_mpoly = mesh->mpoly;
+
+  uint numNewVerts = 0;
+  uint numNewEdges = 0;
+  uint numNewLoops = 0;
+  uint numNewPolys = 0;
+
+#define MOD_SOLIDIFY_EMPTY_TAG ((uint)-1)
+
+  /* Calculate only face normals. */
+  poly_nors = MEM_malloc_arrayN(numPolys, sizeof(*poly_nors), __func__);
+  BKE_mesh_calc_normals_poly(orig_mvert,
+                             NULL,
+                             (int)numVerts,
+                             orig_mloop,
+                             orig_mpoly,
+                             (int)numLoops,
+                             (int)numPolys,
+                             poly_nors,
+                             true);
+
+  NewFaceRef *face_sides_arr = MEM_malloc_arrayN(
+      numPolys * 2, sizeof(*face_sides_arr), "face_sides_arr in solidify");
+  bool *null_faces =
+      (smd->nonmanifold_offset_mode == MOD_SOLIDIFY_NONMANIFOLD_OFFSET_MODE_CONSTRAINTS) ?
+          MEM_calloc_arrayN(numPolys, sizeof(*null_faces), "null_faces in solidify") :
+          NULL;
+  uint largest_ngon = 3;
+  /* Calculate face to #NewFaceRef map. */
+  {
+    mp = orig_mpoly;
+    for (uint i = 0; i < numPolys; i++, mp++) {
+      /* Make normals for faces without area (should really be avoided though). */
+      if (len_squared_v3(poly_nors[i]) < 0.5f) {
+        MEdge *e = orig_medge + orig_mloop[mp->loopstart].e;
+        float edgedir[3];
+        sub_v3_v3v3(edgedir, orig_mvert[e->v2].co, orig_mvert[e->v1].co);
+        if (fabsf(edgedir[2]) < fabsf(edgedir[1])) {
+          poly_nors[i][2] = 1.0f;
+        }
+        else {
+          poly_nors[i][1] = 1.0f;
+        }
+        if (null_faces) {
+          null_faces[i] = true;
+        }
+      }
+
+      NewEdgeRef **link_edges = MEM_calloc_arrayN(
+          (uint)mp->totloop, sizeof(*link_edges), "NewFaceRef::link_edges in solidify");
+      face_sides_arr[i * 2] = (NewFaceRef){
+          .face = mp, .index = i, .reversed = false, .link_edges = link_edges};
+      link_edges = MEM_calloc_arrayN(
+          (uint)mp->totloop, sizeof(*link_edges), "NewFaceRef::link_edges in solidify");
+      face_sides_arr[i * 2 + 1] = (NewFaceRef){
+          .face = mp, .index = i, .reversed = true, .link_edges = link_edges};
+      if (mp->totloop > largest_ngon) {
+        largest_ngon = (uint)mp->totloop;
+      }
+      if (do_shell) {
+        numNewPolys += 2;
+        numNewLoops += (uint)mp->totloop * 2;
+      }
+    }
+  }
+
+  uint *edge_adj_faces_len = MEM_calloc_arrayN(
+      numEdges, sizeof(*edge_adj_faces_len), "edge_adj_faces_len in solidify");
+  /* Count for each edge how many faces it has adjacent. */
+  {
+    mp = orig_mpoly;
+    for (uint i = 0; i < numPolys; i++, mp++) {
+      ml = orig_mloop + mp->loopstart;
+      for (uint j = 0; j < mp->totloop; j++, ml++) {
+        edge_adj_faces_len[ml->e]++;
+      }
+    }
+  }
+
+  /* Original edge to #NewEdgeRef map. */
+  NewEdgeRef ***orig_edge_data_arr = MEM_calloc_arrayN(
+      numEdges, sizeof(*orig_edge_data_arr), "orig_edge_data_arr in solidify");
+  /* Original edge length cache. */
+  float *orig_edge_lengths = MEM_calloc_arrayN(
+      numEdges, sizeof(*orig_edge_lengths), "orig_edge_lengths in solidify");
+  /* Edge groups for every original vert. */
+  EdgeGroup **orig_vert_groups_arr = MEM_calloc_arrayN(
+      numVerts, sizeof(*orig_vert_groups_arr), "orig_vert_groups_arr in solidify");
+  /* Duplicate verts map. */
+  uint *vm = MEM_malloc_arrayN(numVerts, sizeof(*vm), "orig_vert_map in solidify");
+  for (uint i = 0; i < numVerts; i++) {
+    vm[i] = i;
+  }
+
+  uint edge_index = 0;
+  uint loop_index = 0;
+  uint poly_index = 0;
+
+  bool has_singularities = false;
+
+  /* Vert edge adjacent map. */
+  uint *vert_adj_edges_len = MEM_calloc_arrayN(
+      numVerts, sizeof(*vert_adj_edges_len), "vert_adj_edges_len in solidify");
+  OldVertEdgeRef **vert_adj_edges = MEM_calloc_arrayN(
+      numVerts, sizeof(*vert_adj_edges), "vert_adj_edges in solidify");
+
+  /* Create edge to #NewEdgeRef map. */
+  {
+    OldEdgeFaceRef **edge_adj_faces = MEM_calloc_arrayN(
+        numEdges, sizeof(*edge_adj_faces), "edge_adj_faces in solidify");
+
+    /* Create link_faces for edges. */
+    {
+      mp = orig_mpoly;
+      for (uint i = 0; i < numPolys; i++, mp++) {
+        ml = orig_mloop + mp->loopstart;
+        for (uint j = 0; j < mp->totloop; j++, ml++) {
+          const uint edge = ml->e;
+          const bool reversed = orig_medge[edge].v2 != ml->v;
+          OldEdgeFaceRef *old_face_edge_ref = edge_adj_faces[edge];
+          if (old_face_edge_ref == NULL) {
+            const uint len = edge_adj_faces_len[edge];
+            BLI_assert(len > 0);
+            uint *adj_faces = MEM_malloc_arrayN(
+                len, sizeof(*adj_faces), "OldEdgeFaceRef::faces in solidify");
+            bool *adj_faces_loops_reversed = MEM_malloc_arrayN(
+                len, sizeof(*adj_faces_loops_reversed), "OldEdgeFaceRef::reversed in solidify");
+            adj_faces[0] = i;
+            for (uint k = 1; k < len; k++) {
+              adj_faces[k] = MOD_SOLIDIFY_EMPTY_TAG;
+            }
+            adj_faces_loops_reversed[0] = reversed;
+            OldEdgeFaceRef *ref = MEM_mallocN(sizeof(*ref), "OldEdgeFaceRef in solidify");
+            *ref = (OldEdgeFaceRef){adj_faces, len, adj_faces_loops_reversed, 1};
+            edge_adj_faces[edge] = ref;
+          }
+          else {
+            for (uint k = 1; k < old_face_edge_ref->faces_len; k++) {
+              if (old_face_edge_ref->faces[k] == MOD_SOLIDIFY_EMPTY_TAG) {
+                old_face_edge_ref->faces[k] = i;
+                old_face_edge_ref->faces_reversed[k] = reversed;
+                break;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    float edgedir[3] = {0, 0, 0};
+
+    /* Calculate edge lengths and len vert_adj edges. */
+    {
+      bool *face_singularity = MEM_calloc_arrayN(
+          numPolys, sizeof(*face_singularity), "face_sides_arr in solidify");
+      ed = orig_medge;
+      for (uint i = 0; i < numEdges; i++, ed++) {
+        if (edge_adj_faces_len[i] > 0) {
+          const uint v1 = vm[ed->v1];
+          const uint v2 = vm[ed->v2];
+          if (v1 != v2) {
+            sub_v3_v3v3(edgedir, orig_mvert[v2].co, orig_mvert[v1].co);
+            orig_edge_lengths[i] = len_squared_v3(edgedir);
+          }
+          if (v1 == v2 || orig_edge_lengths[i] <= FLT_EPSILON) {
+            if (v2 > v1) {
+              for (uint j = v2; j < numVerts; j++) {
+                if (vm[j] == v2) {
+                  vm[j] = v1;
+                  vert_adj_edges_len[v1] += vert_adj_edges_len[j];
+                  vert_adj_edges_len[j] = 0;
+                }
+              }
+            }
+            else if (v2 < v1) {
+              for (uint j = v1; j < numVerts; j++) {
+                if (vm[j] == v1) {
+                  vm[j] = v2;
+                  vert_adj_edges_len[v2] += vert_adj_edges_len[j];
+                  vert_adj_edges_len[j] = 0;
+                }
+              }
+            }
+            if (do_shell) {
+              numNewLoops -= edge_adj_faces_len[i] * 2;
+            }
+            if (v1 == v2) {
+              /* Remove polys. */
+              for (uint j = 0; j < edge_adj_faces[i]->faces_len; j++) {
+                const uint face = edge_adj_faces[i]->faces[j];
+                if (!face_singularity[face]) {
+                  bool is_singularity = true;
+                  for (uint k = 0; k < orig_mpoly[face].totloop; k++) {
+                    if (vm[orig_mloop[((uint)orig_mpoly[face].loopstart) + k].v] != v1) {
+                      is_singularity = false;
+                      break;
+                    }
+                  }
+                  if (is_singularity) {
+                    face_singularity[face] = true;
+                    if (do_shell) {
+                      numNewPolys -= 2;
+                    }
+                  }
+                }
+              }
+            }
+            edge_adj_faces_len[i] = 0;
+            MEM_freeN(edge_adj_faces[i]->faces);
+            MEM_freeN(edge_adj_faces[i]->faces_reversed);
+            MEM_freeN(edge_adj_faces[i]);
+            edge_adj_faces[i] = NULL;
+          }
+          else if (edge_adj_faces_len[i] > 0) {
+            orig_edge_lengths[i] = sqrtf(orig_edge_lengths[i]);
+            vert_adj_edges_len[v1]++;
+            vert_adj_edges_len[v2]++;
+          }
+        }
+      }
+      MEM_freeN(face_singularity);
+    }
+
+    /* Create vert_adj_edges for verts. */
+    {
+      ed = orig_medge;
+      for (uint i = 0; i < numEdges; i++, ed++) {
+        if (edge_adj_faces_len[i] > 0) {
+          const uint vs[2] = {vm[ed->v1], vm[ed->v2]};
+          uint invalid_edge_index = 0;
+          bool invalid_edge_reversed = false;
+          for (uint j = 0; j < 2; j++) {
+            const uint vert = vs[j];
+            const uint len = vert_adj_edges_len[vert];
+            if (len > 0) {
+              OldVertEdgeRef *old_edge_vert_ref = vert_adj_edges[vert];
+              if (old_edge_vert_ref == NULL) {
+                uint *adj_edges = MEM_calloc_arrayN(
+                    len, sizeof(*adj_edges), "OldVertEdgeRef::edges in solidify");
+                adj_edges[0] = i;
+                for (uint k = 1; k < len; k++) {
+                  adj_edges[k] = MOD_SOLIDIFY_EMPTY_TAG;
+                }
+                OldVertEdgeRef *ref = MEM_mallocN(sizeof(*ref), "OldVertEdgeRef in solidify");
+                *ref = (OldVertEdgeRef){adj_edges, 1};
+                vert_adj_edges[vert] = ref;
+              }
+              else {
+                const uint *f = old_edge_vert_ref->edges;
+                for (uint k = 0; k < len && k <= old_edge_vert_ref->edges_len; k++, f++) {
+                  const uint edge = old_edge_vert_ref->edges[k];
+                  if (edge == MOD_SOLIDIFY_EMPTY_TAG || k == old_edge_vert_ref->edges_len) {
+                    old_edge_vert_ref->edges[k] = i;
+                    old_edge_vert_ref->edges_len++;
+                    break;
+                  }
+                  else if (vm[orig_medge[edge].v1] == vs[1 - j]) {
+                    invalid_edge_index = edge + 1;
+                    invalid_edge_reversed = (j == 0);
+                    break;
+                  }
+                  else if (vm[orig_medge[edge].v2] == vs[1 - j]) {
+                    invalid_edge_index = edge + 1;
+                    invalid_edge_reversed = (j == 1);
+                    break;
+                  }
+                }
+                if (invalid_edge_index) {
+                  /* Should never actually be executed. */
+                  if (j == 1) {
+                    vert_adj_edges[vs[0]]->edges_len--;
+                  }
+                  break;
+                }
+              }
+            }
+          }
+          if (invalid_edge_index) {
+            const uint tmp = invalid_edge_index - 1;
+            invalid_edge_index = i;
+            i = tmp;
+            OldEdgeFaceRef *i_adj_faces = edge_adj_faces[i];
+            OldEdgeFaceRef *invalid_adj_faces = edge_adj_faces[invalid_edge_index];
+            uint j = 0;
+            for (uint k = 0; k < i_adj_faces->faces_len; k++) {
+              for (uint l = 0; l < invalid_adj_faces->faces_len; l++) {
+                if (i_adj_faces->faces[k] == invalid_adj_faces->faces[l] &&
+                    i_adj_faces->faces[k] != MOD_SOLIDIFY_EMPTY_TAG) {
+                  i_adj_faces->faces[k] = MOD_SOLIDIFY_EMPTY_TAG;
+                  invalid_adj_faces->faces[l] = MOD_SOLIDIFY_EMPTY_TAG;
+                  j++;
+                }
+              }
+            }
+            if (do_shell) {
+              numNewPolys -= 2 * j;
+              numNewLoops -= 4 * j;
+            }
+            const uint len = i_adj_faces->faces_len + invalid_adj_faces->faces_len - 2 * j;
+            uint *adj_faces = MEM_malloc_arrayN(
+                len, sizeof(*adj_faces), "OldEdgeFaceRef::faces in solidify");
+            bool *adj_faces_loops_reversed = MEM_malloc_arrayN(
+                len, sizeof(*adj_faces_loops_reversed), "OldEdgeFaceRef::reversed in solidify");
+            j = 0;
+            for (uint k = 0; k < i_adj_faces->faces_len; k++) {
+              if (i_adj_faces->faces[k] != MOD_SOLIDIFY_EMPTY_TAG) {
+                adj_faces[j] = i_adj_faces->faces[k];
+                adj_faces_loops_reversed[j++] = i_adj_faces->faces_reversed[k];
+              }
+            }
+            for (uint k = 0; k < invalid_adj_faces->faces_len; k++) {
+              if (invalid_adj_faces->faces[k] != MOD_SOLIDIFY_EMPTY_TAG) {
+                adj_faces[j] = invalid_adj_faces->faces[k];
+                adj_faces_loops_reversed[j++] = (invalid_edge_reversed !=
+                                                 invalid_adj_faces->faces_reversed[k]);
+              }
+            }
+            BLI_assert(j == len);
+            edge_adj_faces_len[invalid_edge_index] = 0;
+            edge_adj_faces_len[i] = len;
+            MEM_freeN(i_adj_faces->faces);
+            MEM_freeN(i_adj_faces->faces_reversed);
+            i_adj_faces->faces_len = len;
+            i_adj_faces->faces = adj_faces;
+            i_adj_faces->faces_reversed = adj_faces_loops_reversed;
+            i_adj_faces->used += invalid_adj_faces->used;
+            MEM_freeN(invalid_adj_faces->faces);
+            MEM_freeN(invalid_adj_faces->faces_reversed);
+            MEM_freeN(invalid_adj_faces);
+            edge_adj_faces[invalid_edge_index] = i_adj_faces;
+            /* Reset counter to continue. */
+            i = invalid_edge_index;
+          }
+        }
+      }
+    }
+
+    /* Filter duplicate polys. */
+    {
+      ed = orig_medge;
+      for (uint i = 0; i < numEdges; i++, ed++) {
+        if (edge_adj_faces_len[i] > 0) {
+          const OldEdgeFaceRef *adj_faces = edge_adj_faces[i];
+          uint adj_len = adj_faces->faces_len;
+          if (adj_len > 1) {
+            /* For each face pair check if they have equal verts. */
+            for (uint j = 0; j < adj_len; j++) {
+              const uint face = adj_faces->faces[j];
+              const int j_loopstart = orig_mpoly[face].loopstart;
+              const int totloop = orig_mpoly[face].totloop;
+              const uint j_first_v = vm[orig_mloop[j_loopstart].v];
+              for (uint k = j + 1; k < adj_len; k++) {
+                if (orig_mpoly[adj_faces->faces[k]].totloop != totloop) {
+                  continue;
+                }
+                /* Find first face first loop vert in second face loops. */
+                const int k_loopstart = orig_mpoly[adj_faces->faces[k]].loopstart;
+                int l;
+                ml = orig_mloop + k_loopstart;
+                for (l = 0; l < totloop && vm[ml->v] != j_first_v; l++, ml++) {
+                  /* Pass. */
+                }
+                if (l == totloop) {
+                  continue;
+                }
+                /* Check if all following loops have equal verts. */
+                const bool reversed = adj_faces->faces_reversed[j] != adj_faces->faces_reversed[k];
+                const int count_dir = reversed ? -1 : 1;
+                bool has_diff = false;
+                ml = orig_mloop + j_loopstart;
+                for (int m = 0, n = l + totloop; m < totloop && !has_diff;
+                     m++, n += count_dir, ml++) {
+                  has_diff = has_diff || vm[ml->v] != vm[orig_mloop[k_loopstart + n % totloop].v];
+                }
+                /* If the faces are equal, discard one (j). */
+                if (!has_diff) {
+                  ml = orig_mloop + j_loopstart;
+                  uint del_loops = 0;
+                  for (uint m = 0; m < totloop; m++, ml++) {
+                    const uint e = ml->e;
+                    uint face_index = j;
+                    uint *e_adj_faces_faces = edge_adj_faces[e]->faces;
+                    bool *e_adj_faces_reversed = edge_adj_faces[e]->faces_reversed;
+                    const uint faces_len = edge_adj_faces[e]->faces_len;
+                    if (e != i) {
+                      /* Find index of e in #adj_faces. */
+                      for (face_index = 0;
+                           face_index < faces_len && e_adj_faces_faces[face_index] != face;
+                           face_index++) {
+                        /* Pass. */
+                      }
+                      /* If not found. */
+                      if (face_index == faces_len) {
+                        continue;
+                      }
+                    }
+                    else {
+                      adj_len--;
+                    }
+                    memmove(e_adj_faces_faces + face_index,
+                            e_adj_faces_faces + face_index + 1,
+                            (faces_len - face_index - 1) * sizeof(*e_adj_faces_faces));
+                    memmove(e_adj_faces_reversed + face_index,
+                            e_adj_faces_reversed + face_index + 1,
+                            (faces_len - face_index - 1) * sizeof(*e_adj_faces_reversed));
+                    edge_adj_faces[e]->faces_len--;
+                    if (edge_adj_faces_len[e] > 0) {
+                      edge_adj_faces_len[e]--;
+                      if (edge_adj_faces_len[e] == 0) {
+                        edge_adj_faces[e]->used--;
+                        edge_adj_faces[e] = NULL;
+                      }
+                    }
+                    del_loops++;
+                  }
+                  if (do_shell) {
+                    numNewPolys -= 2;
+                    numNewLoops -= 2 * (uint)del_loops;
+                  }
+                  break;
+                }
+              }
+            }
+          }
+        }
+      }
+    }
+
+    /* Create #NewEdgeRef array. */
+    {
+      ed = orig_medge;
+      for (uint i = 0; i < numEdges; i++, ed++) {
+        const uint v1 = vm[ed->v1];
+        const uint v2 = vm[ed->v2];
+        if (edge_adj_faces_len[i] > 0) {
+          sub_v3_v3v3(edgedir, orig_mvert[v2].co, orig_mvert[v1].co);
+          mul_v3_fl(edgedir, 1.0f / orig_edge_lengths[i]);
+
+          OldEdgeFaceRef *adj_faces = edge_adj_faces[i];
+          const uint adj_len = adj_faces->faces_len;
+          const uint *adj_faces_faces = adj_faces->faces;
+          const bool *adj_faces_reversed = adj_faces->faces_reversed;
+          uint new_edges_len = 0;
+          FaceKeyPair *sorted_faces = MEM_malloc_arrayN(
+              adj_len, sizeof(*sorted_faces), "sorted_faces in solidify");
+          if (adj_len > 1) {
+            new_edges_len = adj_len;
+            /* Get keys for sorting. */
+            float ref_nor[3] = {0, 0, 0};
+            float nor[3];
+            for (uint j = 0; j < adj_len; j++) {
+              const bool reverse = adj_faces_reversed[j];
+              const uint face_i = adj_faces_faces[j];
+              if (reverse) {
+                negate_v3_v3(nor, poly_nors[face_i]);
+              }
+              else {
+                copy_v3_v3(nor, poly_nors[face_i]);
+              }
+              float d = 1;
+              if (orig_mpoly[face_i].totloop > 3) {
+                d = project_v3_v3(nor, edgedir);
+                if (LIKELY(d != 0)) {
+                  d = normalize_v3(nor);
+                }
+                else {
+                  d = 1;
+                }
+              }
+              if (UNLIKELY(d == 0.0f)) {
+                sorted_faces[j].angle = 0.0f;
+              }
+              else if (j == 0) {
+                copy_v3_v3(ref_nor, nor);
+                sorted_faces[j].angle = 0.0f;
+              }
+              else {
+                float angle = angle_signed_on_axis_normalized_v3v3_v3(nor, ref_nor, edgedir);
+                sorted_faces[j].angle = -angle;
+              }
+              sorted_faces[j].face = face_sides_arr + adj_faces_faces[j] * 2 +
+                                     (adj_faces_reversed[j] ? 1 : 0);
+            }
+            /* Sort faces by order around the edge (keep order in faces,
+             * reversed and face_angles the same). */
+            qsort(sorted_faces, adj_len, sizeof(*sorted_faces), comp_float_int_pair);
+          }
+          else {
+            new_edges_len = 2;
+            sorted_faces[0].face = face_sides_arr + adj_faces_faces[0] * 2 +
+                                   (adj_faces_reversed[0] ? 1 : 0);
+            if (do_rim) {
+              /* Only add the loops parallel to the edge for now. */
+              numNewLoops += 2;
+              numNewPolys++;
+            }
+          }
+
+          /* Create a list of new edges and fill it. */
+          NewEdgeRef **new_edges = MEM_malloc_arrayN(
+              new_edges_len + 1, sizeof(*new_edges), "new_edges in solidify");
+          new_edges[new_edges_len] = NULL;
+          NewFaceRef *faces[2];
+          for (uint j = 0; j < new_edges_len; j++) {
+            float angle;
+            if (adj_len > 1) {
+              const uint next_j = j + 1 == adj_len ? 0 : j + 1;
+              faces[0] = sorted_faces[j].face;
+              faces[1] = sorted_faces[next_j].face->reversed ? sorted_faces[next_j].face - 1 :
+                                                               sorted_faces[next_j].face + 1;
+              angle = sorted_faces[next_j].angle - sorted_faces[j].angle;
+              if (angle < 0) {
+                angle += 2 * M_PI;
+              }
+            }
+            else {
+              faces[0] = sorted_faces[0].face->reversed ? sorted_faces[0].face - j :
+                                                          sorted_faces[0].face + j;
+              faces[1] = NULL;
+              angle = 0;
+            }
+            NewEdgeRef *edge_data = MEM_mallocN(sizeof(*edge_data), "edge_data in solidify");
+            uint edge_data_edge_index = MOD_SOLIDIFY_EMPTY_TAG;
+            if (do_shell || (adj_len == 1 && do_rim)) {
+              edge_data_edge_index = 0;
+            }
+            *edge_data = (NewEdgeRef){.old_edge = i,
+                                      .faces = {faces[0], faces[1]},
+                                      .link_edge_groups = {NULL, NULL},
+                                      .angle = angle,
+                                      .new_edge = edge_data_edge_index};
+            new_edges[j] = edge_data;
+            for (uint k = 0; k < 2; k++) {
+              if (faces[k] != NULL) {
+                ml = orig_mloop + faces[k]->face->loopstart;
+                for (int l = 0; l < faces[k]->face->totloop; l++, ml++) {
+                  if (edge_adj_faces[ml->e] == edge_adj_faces[i]) {
+                    if (ml->e != i && orig_edge_data_arr[ml->e] == NULL) {
+                      orig_edge_data_arr[ml->e] = new_edges;
+                    }
+                    faces[k]->link_edges[l] = edge_data;
+                    break;
+                  }
+                }
+              }
+            }
+          }
+          MEM_freeN(sorted_faces);
+          orig_edge_data_arr[i] = new_edges;
+          if (do_shell || (adj_len == 1 && do_rim)) {
+            numNewEdges += new_edges_len;
+          }
+        }
+      }
+    }
+
+    for (uint i = 0; i < numEdges; i++) {
+      if (edge_adj_faces[i]) {
+        if (edge_adj_faces[i]->used > 1) {
+          edge_adj_faces[i]->used--;
+        }
+        else {
+          MEM_freeN(edge_adj_faces[i]->faces);
+          MEM_freeN(edge_adj_faces[i]->faces_reversed);
+          MEM_freeN(edge_adj_faces[i]);
+        }
+      }
+    }
+    MEM_freeN(edge_adj_faces);
+  }
+
+  MEM_freeN(vert_adj_edges_len);
+
+  /* Create sorted edge groups for every vert. */
+  {
+    OldVertEdgeRef **adj_edges_ptr = vert_adj_edges;
+    for (uint i = 0; i < numVerts; i++, adj_edges_ptr++) {
+      if (*adj_edges_ptr != NULL && (*adj_edges_ptr)->edges_len >= 2) {
+        EdgeGroup *edge_groups;
+
+        int eg_index = -1;
+        bool contains_long_groups = false;
+        uint topo_groups = 0;
+
+        /* Initial sorted creation. */
+        {
+          const uint *adj_edges = (*adj_edges_ptr)->edges;
+          const uint tot_adj_edges = (*adj_edges_ptr)->edges_len;
+
+          uint unassigned_edges_len = 0;
+          for (uint j = 0; j < tot_adj_edges; j++) {
+            NewEdgeRef **new_edges = orig_edge_data_arr[adj_edges[j]];
+            /* TODO check where the null pointer come from,
+             * because there should not be any... */
+            if (new_edges) {
+              while (*new_edges) {
+                unassigned_edges_len++;
+                new_edges++;
+              }
+            }
+          }
+          NewEdgeRef **unassigned_edges = MEM_malloc_arrayN(
+              unassigned_edges_len, sizeof(*unassigned_edges), "unassigned_edges in solidify");
+          for (uint j = 0, k = 0; j < tot_adj_edges; j++) {
+            NewEdgeRef **new_edges = orig_edge_data_arr[adj_edges[j]];
+            if (new_edges) {
+              while (*new_edges) {
+                unassigned_edges[k++] = *new_edges;
+                new_edges++;
+              }
+            }
+          }
+
+          edge_groups = MEM_calloc_arrayN(
+              (unassigned_edges_len / 2) + 1, sizeof(*edge_groups), "edge_groups in solidify");
+
+          uint assigned_edges_len = 0;
+          NewEdgeRef *found_edge = NULL;
+          uint found_edge_index = 0;
+          bool insert_at_start = false;
+          uint eg_capacity = 5;
+          NewFaceRef *eg_track_faces[2] = {NULL, NULL};
+          NewFaceRef *last_open_edge_track = NULL;
+          NewEdgeRef *edge = NULL;
+
+          while (assigned_edges_len < unassigned_edges_len) {
+            found_edge = NULL;
+            insert_at_start = false;
+            if (eg_index >= 0 && edge_groups[eg_index].edges_len == 0) {
+              uint j = 0;
+              edge = NULL;
+              while (!edge && j < unassigned_edges_len) {
+                edge = unassigned_edges[j++];
+                if (edge && last_open_edge_track &&
+                    (edge->faces[0] != last_open_edge_track || edge->faces[1] != NULL)) {
+                  edge = NULL;
+                }
+              }
+              if (!edge && last_open_edge_track) {
+                topo_groups++;
+                last_open_edge_track = NULL;
+                edge_groups[eg_index].topo_group++;
+                j = 0;
+                while (!edge && j < unassigned_edges_len) {
+                  edge = unassigned_edges[j++];
+                }
+              }
+              else if (!last_open_edge_track && eg_index > 0) {
+                topo_groups++;
+                edge_groups[eg_index].topo_group++;
+              }
+              BLI_assert(edge != NULL);
+              found_edge_index = j - 1;
+              found_edge = edge;
+              if (!last_open_edge_track && vm[orig_medge[edge->old_edge].v1] == i) {
+                eg_track_faces[0] = edge->faces[0];
+                eg_track_faces[1] = edge->faces[1];
+                if (edge->faces[1] == NULL) {
+                  last_open_edge_track = edge->faces[0]->reversed ? edge->faces[0] - 1 :
+                                                                    edge->faces[0] + 1;
+                }
+              }
+              else {
+                eg_track_faces[0] = edge->faces[1];
+                eg_track_faces[1] = edge->faces[0];
+              }
+            }
+            else if (eg_index >= 0) {
+              NewEdgeRef **edge_ptr = unassigned_edges;
+              for (found_edge_index = 0; found_edge_index < unassigned_edges_len;
+                   found_edge_index++, edge_ptr++) {
+                if (*edge_ptr) {
+                  edge = *edge_ptr;
+                  if (edge->faces[0] == eg_track_faces[1]) {
+                    insert_at_start = false;
+                    eg_track_faces[1] = edge->faces[1];
+                    found_edge = edge;
+                    if (edge->faces[1] == NULL) {
+                      edge_groups[eg_index].is_orig_closed = false;
+                      last_open_edge_track = edge->faces[0]->reversed ? edge->faces[0] - 1 :
+                                                                        edge->faces[0] + 1;
+                    }
+                    break;
+                  }
+                  else if (edge->faces[0] == eg_track_faces[0]) {
+                    insert_at_start = true;
+                    eg_track_faces[0] = edge->faces[1];
+                    found_edge = edge;
+                    if (edge->faces[1] == NULL) {
+                      edge_groups[eg_index].is_orig_closed = false;
+                    }
+                    break;
+                  }
+                  else if (edge->faces[1] != NULL) {
+                    if (edge->faces[1] == eg_track_faces[1]) {
+                      insert_at_start = false;
+                      eg_track_faces[1] = edge->faces[0];
+                      found_edge = edge;
+                      break;
+                    }
+                    else if (edge->faces[1] == eg_track_faces[0]) {
+                      insert_at_start = true;
+                      eg_track_faces[0] = edge->faces[0];
+                      found_edge = edge;
+                      break;
+                    }
+                  }
+                }
+              }
+            }
+            if (found_edge) {
+              unassigned_edges[found_edge_index] = NULL;
+              assigned_edges_len++;
+              const uint needed_capacity = edge_groups[eg_index].edges_len + 1;
+              if (needed_capacity > eg_capacity) {
+                eg_capacity = needed_capacity + 1;
+                NewEdgeRef **new_eg = MEM_calloc_arrayN(
+                    eg_capacity, sizeof(*new_eg), "edge_group realloc in solidify");
+                if (insert_at_start) {
+                  memcpy(new_eg + 1,
+                         edge_groups[eg_index].edges,
+                         edge_groups[eg_index].edges_len * sizeof(*new_eg));
+                }
+                else {
+                  memcpy(new_eg,
+                         edge_groups[eg_index].edges,
+                         edge_groups[eg_index].edges_len * sizeof(*new_eg));
+                }
+                MEM_freeN(edge_groups[eg_index].edges);
+                edge_groups[eg_index].edges = new_eg;
+              }
+              else if (insert_at_start) {
+                memmove(edge_groups[eg_index].edges + 1,
+                        edge_groups[eg_index].edges,
+                        edge_groups[eg_index].edges_len * sizeof(*edge_groups[eg_index].edges));
+              }
+              edge_groups[eg_index].edges[insert_at_start ? 0 : edge_groups[eg_index].edges_len] =
+                  found_edge;
+              edge_groups[eg_index].edges_len++;
+              if (edge_groups[eg_index].edges[edge_groups[eg_index].edges_len - 1]->faces[1] !=
+                  NULL) {
+                last_open_edge_track = NULL;
+              }
+              if (edge_groups[eg_index].edges_len > 3) {
+                contains_long_groups = true;
+              }
+            }
+            else {
+              eg_index++;
+              BLI_assert(eg_index < (unassigned_edges_len / 2));
+              eg_capacity = 5;
+              NewEdgeRef **edges = MEM_calloc_arrayN(
+                  eg_capacity, sizeof(*edges), "edge_group in solidify");
+              edge_groups[eg_index] = (EdgeGroup){
+                  .valid = true,
+                  .edges = edges,
+                  .edges_len = 0,
+                  .open_face_edge = MOD_SOLIDIFY_EMPTY_TAG,
+                  .is_orig_closed = true,
+                  .is_even_split = false,
+                  .split = 0,
+                  .is_singularity = false,
+                  .topo_group = topo_groups,
+                  .co = {0.0f, 0.0f, 0.0f},
+                  .no = {0.0f, 0.0f, 0.0f},
+                  .new_vert = MOD_SOLIDIFY_EMPTY_TAG,
+              };
+              eg_track_faces[0] = NULL;
+              eg_track_faces[1] = NULL;
+            }
+          }
+          /* #eg_index is the number of groups from here on. */
+          eg_index++;
+          /* #topo_groups is the number of topo groups from here on. */
+          topo_groups++;
+          MEM_freeN(unassigned_edges);
+        }
+
+        /* Split of long self intersection groups */
+        {
+          uint splits = 0;
+          if (contains_long_groups) {
+            uint add_index = 0;
+            for (uint j = 0; j < eg_index; j++) {
+              const uint edges_len = edge_groups[j + add_index].edges_len;
+              if (edges_len > 3) {
+                bool has_doubles = false;
+                bool *doubles = MEM_calloc_arrayN(
+                    edges_len, sizeof(*doubles), "doubles in solidify");
+                EdgeGroup g = edge_groups[j + add_index];
+                for (uint k = 0; k < edges_len; k++) {
+                  for (uint l = k + 1; l < edges_len; l++) {
+                    if (g.edges[k]->old_edge == g.edges[l]->old_edge) {
+                      doubles[k] = true;
+                      doubles[l] = true;
+                      has_doubles = true;
+                    }
+                  }
+                }
+                if (has_doubles) {
+                  const uint prior_splits = splits;
+                  const uint prior_index = add_index;
+                  int unique_start = -1;
+                  int first_unique_end = -1;
+                  int last_split = -1;
+                  int first_split = -1;
+                  bool first_even_split = false;
+                  uint real_k = 0;
+                  while (real_k < edges_len ||
+                         (g.is_orig_closed &&
+                          (real_k <=
+                               (first_unique_end == -1 ? 0 : first_unique_end) + (int)edges_len ||
+                           first_split != last_split))) {
+                    const uint k = real_k % edges_len;
+                    if (!doubles[k]) {
+                      if (first_unique_end != -1 && unique_start == -1) {
+                        unique_start = (int)real_k;
+                      }
+                    }
+                    else if (first_unique_end == -1) {
+                      first_unique_end = (int)k;
+                    }
+                    else if (unique_start != -1) {
+                      const uint split = (((uint)unique_start + real_k + 1) / 2) % edges_len;
+                      const bool is_even_split = (((uint)unique_start + real_k) & 1);
+                      if (last_split != -1) {
+                        /* Override g on first split (no insert). */
+                        if (prior_splits != splits) {
+                          memmove(edge_groups + j + add_index + 1,
+                                  edge_groups + j + add_index,
+                                  ((uint)eg_index - j) * sizeof(*edge_groups));
+                          add_index++;
+                        }
+                        if (last_split > split) {
+                          const uint size = (split + edges_len) - (uint)last_split;
+                          NewEdgeRef **edges = MEM_malloc_arrayN(
+                              size, sizeof(*edges), "edge_group split in solidify");
+                          memcpy(edges,
+                                 g.edges + last_split,
+                                 (edges_len - (uint)last_split) * sizeof(*edges));
+                          memcpy(edges + (edges_len - (uint)last_split),
+                                 g.edges,
+                                 split * sizeof(*edges));
+                          edge_groups[j + add_index] = (EdgeGroup){
+                              .valid = true,
+                              .edges = edges,
+                              .edges_len = size,
+                              .open_face_edge = MOD_SOLIDIFY_EMPTY_TAG,
+                              .is_orig_closed = g.is_orig_closed,
+                              .is_even_split = is_even_split,
+                              .split = add_index - prior_index + 1 + (uint)!g.is_orig_closed,
+                              .is_singularity = false,
+                              .topo_group = g.topo_group,
+                              .co = {0.0f, 0.0f, 0.0f},
+                              .no = {0.0f, 0.0f, 0.0f},
+                              .new_vert = MOD_SOLIDIFY_EMPTY_TAG,
+                          };
+                        }
+                        else {
+                          const uint size = split - (uint)last_split;
+                          NewEdgeRef **edges = MEM_malloc_arrayN(
+                              size, sizeof(*edges), "edge_group split in solidify");
+                          memcpy(edges, g.edges + last_split, size * sizeof(*edges));
+                          edge_groups[j + add_index] = (EdgeGroup){
+                              .valid = true,
+                              .edges = edges,
+                              .edges_len = size,
+                              .open_face_edge = MOD_SOLIDIFY_EMPTY_TAG,
+                              .is_orig_closed = g.is_orig_closed,
+                              .is_even_split = is_even_split,
+                              .split = add_index - prior_index + 1 + (uint)!g.is_orig_closed,
+                              .is_singularity = false,
+                              .topo_group = g.topo_group,
+                              .co = {0.0f, 0.0f, 0.0f},
+                              .no = {0.0f, 0.0f, 0.0f},
+                              .new_vert = MOD_SOLIDIFY_EMPTY_TAG,
+                          };
+                        }
+                        splits++;
+                      }
+                      last_split = (int)split;
+                      if (first_split == -1) {
+                        first_split = (int)split;
+                        first_even_split = is_even_split;
+                      }
+                      unique_start = -1;
+                    }
+                    real_k++;
+                  }
+                  if (first_split != -1) {
+                    if (!g.is_orig_closed) {
+                      if (prior_splits != splits) {
+                        memmove(edge_groups + (j + prior_index + 1),
+                                edge_groups + (j + prior_index),
+                                ((uint)eg_index + add_index - (j + prior_index)) *
+                                    sizeof(*edge_groups));
+                        memmove(edge_groups + (j + add_index + 2),
+                                edge_groups + (j + add_index + 1),
+                                ((uint)eg_index - j) * sizeof(*edge_groups));
+                        add_index++;
+                      }
+                      else {
+                        memmove(edge_groups + (j + add_index + 2),
+                                edge_groups + (j + add_index + 1),
+                                ((uint)eg_index - j - 1) * sizeof(*edge_groups));
+                      }
+                      NewEdgeRef **edges = MEM_malloc_arrayN(
+                          (uint)first_split, sizeof(*edges), "edge_group split in solidify");
+                      memcpy(edges, g.edges, (uint)first_split * sizeof(*edges));
+                      edge_groups[j + prior_index] = (EdgeGroup){
+                          .valid = true,
+                          .edges = edges,
+                          .edges_len = (uint)first_split,
+                          .open_face_edge = MOD_SOLIDIFY_EMPTY_TAG,
+                          .is_orig_closed = g.is_orig_closed,
+                          .is_even_split = first_even_split,
+                          .split = 1,
+                          .is_singularity = false,
+                          .topo_group = g.topo_group,
+                          .co = {0.0f, 0.0f, 0.0f},
+                          .no = {0.0f, 0.0f, 0.0f},
+                          .new_vert = MOD_SOLIDIFY_EMPTY_TAG,
+                      };
+                      add_index++;
+                      splits++;
+                      edges = MEM_malloc_arrayN(edges_len - (uint)last_split,
+                                                sizeof(*edges),
+                                                "edge_group split in solidify");
+                      memcpy(edges,
+                             g.edges + last_split,
+                             (edges_len - (uint)last_split) * sizeof(*edges));
+                      edge_groups[j + add_index] = (EdgeGroup){
+                          .valid = true,
+                          .edges = edges,
+                          .edges_len = (edges_len - (uint)last_split),
+                          .open_face_edge = MOD_SOLIDIFY_EMPTY_TAG,
+                          .is_orig_closed = g.is_orig_closed,
+                          .is_even_split = false,
+                          .split = add_index - prior_index + 1,
+                          .is_singularity = false,
+                          .topo_group = g.topo_group,
+                          .co = {0.0f, 0.0f, 0.0f},
+                          .no = {0.0f, 0.0f, 0.0f},
+                          .new_vert = MOD_SOLIDIFY_EMPTY_TAG,
+                      };
+                    }
+                    if (prior_splits != splits) {
+                      MEM_freeN(g.edges);
+                    }
+                  }
+                  if (first_unique_end != -1 && prior_splits == splits) {
+                    has_singularities = true;
+                    edge_groups[j + add_index].is_singularity = true;
+                  }
+                }
+                MEM_freeN(doubles);
+              }
+            }
+          }
+        }
+
+        orig_vert_groups_arr[i] = edge_groups;
+        /* Count new edges, loops, polys and add to link_edge_groups. */
+        {
+          uint new_verts = 0;
+          bool contains_open_splits = false;
+          uint open_edges = 0;
+          uint contains_splits = 0;
+          uint last_added = 0;
+          uint first_added = 0;
+          bool first_set = false;
+          for (EdgeGroup *g = edge_groups; g->valid; g++) {
+            NewEdgeRef **e = g->edges;
+            for (uint j = 0; j < g->edges_len; j++, e++) {
+              const uint flip = (uint)(vm[orig_medge[(*e)->old_edge].v2] == i);
+              BLI_assert(flip || vm[orig_medge[(*e)->old_edge].v1] == i);
+              (*e)->link_edge_groups[flip] = g;
+            }
+            uint added = 0;
+            if (do_shell || (do_rim && !g->is_orig_closed)) {
+              BLI_assert(g->new_vert == MOD_SOLIDIFY_EMPTY_TAG);
+              g->new_vert = numNewVerts++;
+              if (do_rim || (do_shell && g->split)) {
+                new_verts++;
+                contains_splits += (g->split != 0);
+                contains_open_splits |= g->split && !g->is_orig_closed;
+                added = g->split;
+              }
+            }
+            open_edges += (uint)(added < last_added);
+            if (!first_set) {
+              first_set = true;
+              first_added = added;
+            }
+            last_added = added;
+            if (!(g + 1)->valid || g->topo_group != (g + 1)->topo_group) {
+              if (new_verts > 2) {
+                numNewPolys++;
+                numNewEdges += new_verts;
+                open_edges += (uint)(first_added < last_added);
+                open_edges -= (uint)(open_edges && !contains_open_splits);
+                if (do_shell && do_rim) {
+                  numNewLoops += new_verts * 2;
+                }
+                else if (do_shell) {
+                  numNewLoops += new_verts * 2 - open_edges;
+                }
+                else {  // do_rim
+                  numNewLoops += new_verts * 2 + open_edges - contains_splits;
+                }
+              }
+              else if (new_verts == 2) {
+                numNewEdges++;
+                numNewLoops += 2u - (uint)(!(do_rim && do_shell) && contains_open_splits);
+              }
+              new_verts = 0;
+              contains_open_splits = false;
+              contains_splits = 0;
+              open_edges = 0;
+              last_added = 0;
+              first_added = 0;
+              first_set = false;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  /* Free vert_adj_edges memory. */
+  {
+    uint i = 0;
+    for (OldVertEdgeRef **p = vert_adj_edges; i < numVerts; i++, p++) {
+      if (*p) {
+        MEM_freeN((*p)->edges);
+        MEM_freeN(*p);
+      }
+    }
+    MEM_freeN(vert_adj_edges);
+  }
+
+  /* TODO create_regions if fix_intersections. */
+
+  /* General use pointer for #EdgeGroup iteration. */
+  EdgeGroup **gs_ptr;
+
+  /* Calculate EdgeGroup vertex coordinates. */
+  {
+    mv = orig_mvert;
+    gs_ptr = orig_vert_groups_arr;
+    for (uint i = 0; i < numVerts; i++, mv++, gs_ptr++) {
+      if (*gs_ptr) {
+        EdgeGroup *g = *gs_ptr;
+        for (uint j = 0; g->valid; j++, g++) {
+          if (!g->is_singularity) {
+            float *nor = g->no;
+            float move_nor[3] = {0, 0, 0};
+            bool disable_boundary_fix = (smd->nonmanifold_boundary_mode ==
+                                             MOD_SOLIDIFY_NONMANIFOLD_BOUNDARY_MODE_NONE ||
+                                         (g->is_orig_closed || g->split));
+            /* Constraints Method. */
+            if (smd->nonmanifold_offset_mode == MOD_SOLIDIFY_NONMANIFOLD_OFFSET_MODE_CONSTRAINTS) {
+              NewEdgeRef *first_edge = NULL;
+              NewEdgeRef **edge_ptr = g->edges;
+              /* Contains normal and offset [nx, ny, nz, ofs]. */
+              float(*normals_queue)[4] = MEM_malloc_arrayN(
+                  g->edges_len + 1, sizeof(*normals_queue), "normals_queue in solidify");
+              uint queue_index = 0;
+
+              float face_nors[3][3];
+              float nor_ofs[3];
+
+              const bool cycle = (g->is_orig_closed && !g->split) || g->is_even_split;
+              for (uint k = 0; k < g->edges_len; k++, edge_ptr++) {
+                if (!(k & 1) || (!cycle && k == g->edges_len - 1)) {
+                  NewEdgeRef *edge = *edge_ptr;
+                  for (uint l = 0; l < 2; l++) {
+                    NewFaceRef *face = edge->faces[l];
+                    if (face && (first_edge == NULL ||
+                                 (first_edge->faces[0] != face && first_edge->faces[1] != face))) {
+                      if (!null_faces[face->index]) {
+                        mul_v3_v3fl(normals_queue[queue_index],
+                                    poly_nors[face->index],
+                                    face->reversed ? -1 : 1);
+                        normals_queue[queue_index++][3] = face->reversed ? ofs_back : ofs_front;
+                      }
+                      else {
+                        mul_v3_v3fl(face_nors[0], poly_nors[face->index], face->reversed ? -1 : 1);
+                        nor_ofs[0] = face->reversed ? ofs_back : ofs_front;
+                      }
+                    }
+                  }
+                  if ((cycle && k == 0) || (!cycle && k + 3 >= g->edges_len)) {
+                    first_edge = edge;
+                  }
+                }
+              }
+              uint face_nors_len = 0;
+              const float stop_explosion = 1 - fabsf(smd->offset_fac) * 0.05f;
+              while (queue_index > 0) {
+                if (face_nors_len == 0) {
+                  if (queue_index <= 2) {
+                    for (uint k = 0; k < queue_index; k++) {
+                      copy_v3_v3(face_nors[k], normals_queue[k]);
+                      nor_ofs[k] = normals_queue[k][3];
+                    }
+                    face_nors_len = queue_index;
+                    queue_index = 0;
+                  }
+                  else {
+                    /* Find most different two normals. */
+                    float min_p = 2;
+                    uint min_n0 = 0;
+                    uint min_n1 = 0;
+                    for (uint k = 0; k < queue_index; k++) {
+                      for (uint m = k + 1; m < queue_index; m++) {
+                        float p = dot_v3v3(normals_queue[k], normals_queue[m]);
+                        if (p <= min_p + FLT_EPSILON) {
+                          min_p = p;
+                          min_n0 = m;
+                          min_n1 = k;
+                        }
+                      }
+                    }
+                    copy_v3_v3(face_nors[0], normals_queue[min_n0]);
+                    copy_v3_v3(face_nors[1], normals_queue[min_n1]);
+                    nor_ofs[0] = normals_queue[min_n0][3];
+                    nor_ofs[1] = normals_queue[min_n1][3];
+                    face_nors_len = 2;
+                    queue_index--;
+                    memmove(normals_queue + min_n0,
+                            normals_queue + min_n0 + 1,
+                            (queue_index - min_n0) * sizeof(*normals_queue));
+                    queue_index--;
+                    memmove(normals_queue + min_n1,
+                            normals_queue + min_n1 + 1,
+                            (queue_index - min_n1) * sizeof(*normals_queue));
+                    min_p = 1;
+                    min_n1 = 0;
+                    float max_p = -1;
+                    for (uint k = 0; k < queue_index; k++) {
+                      max_p = -1;
+                      for (uint m = 0; m < face_nors_len; m++) {
+                        float p = dot_v3v3(face_nors[m], normals_queue[k]);
+                        if (p > max_p + FLT_EPSILON) {
+                          max_p = p;
+                        }
+                      }
+                      if (max_p <= min_p + FLT_EPSILON) {
+                        min_p = max_p;
+                        min_n1 = k;
+                      }
+                    }
+                    if (min_p < 0.8) {
+                      copy_v3_v3(face_nors[2], normals_queue[min_n1]);
+                      nor_ofs[2] = normals_queue[min_n1][3];
+                      face_nors_len++;
+                      queue_index--;
+                      memmove(normals_queue + min_n1,
+                              normals_queue + min_n1 + 1,
+                              (queue_index - min_n1) * sizeof(*normals_queue));
+                    }
+                  }
+                }
+                else {
+                  uint best = 0;
+                  uint best_group = 0;
+                  float best_p = -1.0f;
+                  for (uint k = 0; k < queue_index; k++) {
+                    for (uint m = 0; m < face_nors_len; m++) {
+                      float p = dot_v3v3(face_nors[m], normals_queue[k]);
+                      if (p > best_p + FLT_EPSILON) {
+                        best_p = p;
+                        best = m;
+                        best_group = k;
+                      }
+                    }
+                  }
+                  add_v3_v3(face_nors[best], normals_queue[best_group]);
+                  normalize_v3(face_nors[best]);
+                  nor_ofs[best] = (nor_ofs[best] + normals_queue[best_group][3]) * 0.5f;
+                  queue_index--;
+                  memmove(normals_queue + best_group,
+                          normals_queue + best_group + 1,
+                          (queue_index - best_group) * sizeof(*normals_queue));
+                }
+              }
+              MEM_freeN(normals_queue);
+              /* When up to 3 constraint normals are found. */
+              float d, q;
+              switch (face_nors_len) {
+                case 0:
+                  mul_v3_v3fl(nor, face_nors[0], nor_ofs[0]);
+                  disable_boundary_fix = true;
+                  break;
+                case 1:
+                  mul_v3_v3fl(nor, face_nors[0], nor_ofs[0]);
+                  disable_boundary_fix = true;
+                  break;
+                case 2:
+                  q = dot_v3v3(face_nors[0], face_nors[1]);
+                  d = 1.0f - q * q;
+                  if (LIKELY(d > FLT_EPSILON) && q < stop_explosion) {
+                    d = 1.0f / d;
+                    mul_v3_fl(face_nors[0], (nor_ofs[0] - nor_ofs[1] * q) * d);
+                    mul_v3_fl(face_nors[1], (nor_ofs[1] - nor_ofs[0] * q) * d);
+                    add_v3_v3v3(nor, face_nors[0], face_nors[1]);
+                  }
+                  else {
+                    mul_v3_fl(face_nors[0], nor_ofs[0] * 0.5f);
+                    mul_v3_fl(face_nors[1], nor_ofs[1] * 0.5f);
+                    add_v3_v3v3(nor, face_nors[0], face_nors[1]);
+                  }
+                  if (!disable_boundary_fix) {
+                    cross_v3_v3v3(move_nor, face_nors[0], face_nors[1]);
+                  }
+                  break;
+                case 3:
+                  q = dot_v3v3(face_nors[0], face_nors[1]);
+                  d = 1.0f - q * q;
+                  float *free_nor = move_nor; /* No need to allocate a new array. */
+                  cross_v3_v3v3(free_nor, face_nors[0], face_nors[1]);
+                  if (LIKELY(d > FLT_EPSILON) && q < stop_explosion) {
+                    d = 1.0f / d;
+                    mul_v3_fl(face_nors[0], (nor_ofs[0] - nor_ofs[1] * q) * d);
+                    mul_v3_fl(face_nors[1], (nor_ofs[1] - nor_ofs[0] * q) * d);
+                    add_v3_v3v3(nor, face_nors[0], face_nors[1]);
+                  }
+                  else {
+                    mul_v3_fl(face_nors[0], nor_ofs[0] * 0.5f);
+                    mul_v3_fl(face_nors[1], nor_ofs[1] * 0.5f);
+                    add_v3_v3v3(nor, face_nors[0], face_nors[1]);
+                  }
+                  mul_v3_fl(face_nors[2], nor_ofs[2]);
+                  d = dot_v3v3(face_nors[2], free_nor);
+                  if (LIKELY(fabsf(d) > FLT_EPSILON)) {
+                    sub_v3_v3v3(face_nors[0], nor, face_nors[2]); /* Override face_nor[0]. */
+                    mul_v3_fl(free_nor, dot_v3v3(face_nors[2], face_nors[0]) / d);
+                    sub_v3_v3(nor, free_nor);
+                  }
+                  disable_boundary_fix = true;
+                  break;
+                default:
+                  BLI_assert(0);
+              }
+            }
+            /* Simple/Even Method. */
+            else {
+              float total_angle = 0;
+              float total_angle_back = 0;
+              NewEdgeRef *first_edge = NULL;
+              NewEdgeRef **edge_ptr = g->edges;
+              float face_nor[3];
+              float nor_back[3] = {0, 0, 0};
+              bool has_back = false;
+              bool has_front = false;
+              bool cycle = (g->is_orig_closed && !g->split) || g->is_even_split;
+              for (uint k = 0; k < g->edges_len; k++, edge_ptr++) {
+                if (!(k & 1) || (!cycle && k == g->edges_len - 1)) {
+                  NewEdgeRef *edge = *edge_ptr;
+                  for (uint l = 0; l < 2; l++) {
+                    NewFaceRef *face = edge->faces[l];
+                    if (face && (first_edge == NULL ||
+                                 (first_edge->faces[0] != face && first_edge->faces[1] != face))) {
+                      float angle = 1.0f;
+                      float ofs = face->reversed ? -max_ff(1.0e-5f, ofs_back) :
+                                                   max_ff(1.0e-5f, ofs_front);
+                      if (smd->nonmanifold_offset_mode ==
+                          MOD_SOLIDIFY_NONMANIFOLD_OFFSET_MODE_EVEN) {
+                        MLoop *ml_next = orig_mloop + face->face->loopstart;
+                        ml = ml_next + (face->face->totloop - 1);
+                        MLoop *ml_prev = ml - 1;
+                        for (int m = 0; m < face->face->totloop && vm[ml->v] != i;
+                             m++, ml_next++) {
+                          ml_prev = ml;
+                          ml = ml_next;
+                        }
+                        angle = angle_v3v3v3(
+                            orig_mvert[vm[ml_prev->v]].co, mv->co, orig_mvert[vm[ml_next->v]].co);
+                        if (face->reversed) {
+                          total_angle_back += angle * ofs * ofs;
+                        }
+                        else {
+                          total_angle += angle * ofs * ofs;
+                        }
+                      }
+                      else {
+                        if (face->reversed) {
+                          total_angle_back++;
+                        }
+                        else {
+                          total_angle++;
+                        }
+                      }
+                      mul_v3_v3fl(face_nor, poly_nors[face->index], angle * ofs);
+                      if (face->reversed) {
+                        add_v3_v3(nor_back, face_nor);
+                        has_back = true;
+                      }
+                      else {
+                        add_v3_v3(nor, face_nor);
+                        has_front = true;
+                      }
+                    }
+                  }
+                  if ((cycle && k == 0) || (!cycle && k + 3 >= g->edges_len)) {
+                    first_edge = edge;
+                  }
+                }
+              }
+
+              /* Set normal length with selected method. */
+              if (smd->nonmanifold_offset_mode == MOD_SOLIDIFY_NONMANIFOLD_OFFSET_MODE_EVEN) {
+                float d = dot_v3v3(nor, nor_back);
+                if (has_front) {
+                  float length = len_squared_v3(nor);
+                  if (LIKELY(length > FLT_EPSILON)) {
+                    mul_v3_fl(nor, total_angle / length);
+                  }
+                }
+                if (has_back) {
+                  float length = len_squared_v3(nor_back);
+                  if (LIKELY(length > FLT_EPSILON)) {
+                    mul_v3_fl(nor_back, total_angle_back / length);
+                  }
+                  if (!has_front) {
+                    copy_v3_v3(nor, nor_back);
+                  }
+                }
+                if (has_front && has_back) {
+                  float nor_length = len_v3(nor);
+                  float nor_back_length = len_v3(nor_back);
+                  float q = dot_v3v3(nor, nor_back);
+                  if (LIKELY(fabsf(q) > FLT_EPSILON)) {
+                    q /= nor_length * nor_back_length;
+                  }
+                  d = 1.0f - q * q;
+                  if (LIKELY(d > FLT_EPSILON)) {
+                    d = 1.0f / d;
+                    if (LIKELY(nor_length > FLT_EPSILON)) {
+                      mul_v3_fl(nor, (1 - nor_back_length * q / nor_length) * d);
+                    }
+                    if (LIKELY(nor_back_length > FLT_EPSILON)) {
+                      mul_v3_fl(nor_back, (1 - nor_length * q / nor_back_length) * d);
+                    }
+                    add_v3_v3(nor, nor_back);
+                  }
+                  else {
+                    mul_v3_fl(nor, 0.5f);
+                    mul_v3_fl(nor_back, 0.5f);
+                    add_v3_v3(nor, nor_back);
+                  }
+                }
+              }
+              else {
+                if (has_front && total_angle > FLT_EPSILON) {
+                  mul_v3_fl(nor, 1.0f / total_angle);
+                }
+                if (has_back && total_angle_back > FLT_EPSILON) {
+                  mul_v3_fl(nor_back, 1.0f / total_angle_back);
+                  add_v3_v3(nor, nor_back);
+                }
+              }
+              /* Set move_nor for boundary fix. */
+              if (!disable_boundary_fix && g->edges_len > 2) {
+                edge_ptr = g->edges + 1;
+                float tmp[3];
+                uint k;
+                for (k = 1; k + 1 < g->edges_len; k++, edge_ptr++) {
+                  MEdge *e = orig_medge + (*edge_ptr)->old_edge;
+                  sub_v3_v3v3(tmp, orig_mvert[vm[e->v1] == i ? e->v2 : e->v1].co, mv->co);
+                  add_v3_v3(move_nor, tmp);
+                }
+                if (k == 1) {
+                  disable_boundary_fix = true;
+                }
+                else {
+                  disable_boundary_fix = normalize_v3(move_nor) == 0.0f;
+                }
+              }
+              else {
+                disable_boundary_fix = true;
+              }
+            }
+            /* Fix boundary verts. */
+            if (!disable_boundary_fix) {
+              /* Constraint normal, nor * constr_nor == 0 after this fix. */
+              float constr_nor[3];
+              MEdge *e0_edge = orig_medge + g->edges[0]->old_edge;
+              MEdge *e1_edge = orig_medge + g->edges[g->edges_len - 1]->old_edge;
+              float e0[3];
+              float e1[3];
+              sub_v3_v3v3(
+                  e0, orig_mvert[vm[e0_edge->v1] == i ? e0_edge->v2 : e0_edge->v1].co, mv->co);
+              sub_v3_v3v3(
+                  e1, orig_mvert[vm[e1_edge->v1] == i ? e1_edge->v2 : e1_edge->v1].co, mv->co);
+              if (smd->nonmanifold_boundary_mode == MOD_SOLIDIFY_NONMANIFOLD_BOUNDARY_MODE_FLAT) {
+                cross_v3_v3v3(constr_nor, e0, e1);
+              }
+              else {
+                float f0[3];
+                float f1[3];
+                if (g->edges[0]->faces[0]->reversed) {
+                  negate_v3_v3(f0, poly_nors[g->edges[0]->faces[0]->index]);
+                }
+                else {
+                  copy_v3_v3(f0, poly_nors[g->edges[0]->faces[0]->index]);
+                }
+                if (g->edges[g->edges_len - 1]->faces[0]->reversed) {
+                  negate_v3_v3(f1, poly_nors[g->edges[g->edges_len - 1]->faces[0]->index]);
+                }
+                else {
+                  copy_v3_v3(f1, poly_nors[g->edges[g->edges_len - 1]->faces[0]->index]);
+                }
+                float n0[3];
+                float n1[3];
+                cross_v3_v3v3(n0, e0, f0);
+                cross_v3_v3v3(n1, f1, e1);
+                normalize_v3(n0);
+                normalize_v3(n1);
+                add_v3_v3v3(constr_nor, n0, n1);
+              }
+              float d = dot_v3v3(constr_nor, move_nor);
+              if (LIKELY(fabsf(d) > FLT_EPSILON)) {
+                mul_v3_fl(move_nor, dot_v3v3(constr_nor, nor) / d);
+                sub_v3_v3(nor, move_nor);
+              }
+            }
+            float scalar_vgroup = 1;
+            /* Use vertex group. */
+            if (dvert) {
+              MDeformVert *dv = &dvert[i];
+              if (defgrp_invert) {
+                scalar_vgroup = 1.0f - defvert_find_weight(dv, defgrp_index);
+              }
+              else {
+                scalar_vgroup = defvert_find_weight(dv, defgrp_index);
+              }
+              scalar_vgroup = offset_fac_vg + (scalar_vgroup * offset_fac_vg_inv);
+            }
+            /* Do clamping. */
+            if (do_clamp) {
+              if (do_angle_clamp) {
+                float min_length = 0;
+                float angle = 0.5f * M_PI;
+                uint k = 0;
+                for (NewEdgeRef **p = g->edges; k < g->edges_len; k++, p++) {
+                  float length = orig_edge_lengths[(*p)->old_edge];
+                  float e_ang = (*p)->angle;
+                  if (e_ang > angle) {
+                    angle = e_ang;
+                  }
+                  if (length < min_length || k == 0) {
+                    min_length = length;
+                  }
+                }
+                float cos_ang = cosf(angle * 0.5f);
+                if (cos_ang > 0) {
+                  float max_off = min_length * 0.5f / cos_ang;
+                  if (max_off < offset * 0.5f) {
+                    scalar_vgroup *= max_off / offset * 2;
+                  }
+                }
+              }
+              else {
+                float min_length = 0;
+                uint k = 0;
+                for (NewEdgeRef **p = g->edges; k < g->edges_len; k++, p++) {
+                  float length = orig_edge_lengths[(*p)->old_edge];
+                  if (length < min_length || k == 0) {
+                    min_length = length;
+                  }
+                }
+                if (min_length < offset) {
+                  scalar_vgroup *= min_length / offset;
+                }
+              }
+            }
+            mul_v3_fl(nor, scalar_vgroup);
+            add_v3_v3v3(g->co, nor, mv->co);
+          }
+          else {
+            copy_v3_v3(g->co, mv->co);
+          }
+        }
+      }
+    }
+  }
+
+  if (null_faces) {
+    MEM_freeN(null_faces);
+  }
+
+  /* TODO create vertdata for intersection fixes (intersection fixing per topology region). */
+
+  /* Correction for adjacent one sided groups around a vert to
+   * prevent edge duplicates and null polys. */
+  uint(*singularity_edges)[2] = NULL;
+  uint totsingularity = 0;
+  if (has_singularities) {
+    has_singularities = false;
+    uint i = 0;
+    uint singularity_edges_len = 1;
+    singularity_edges = MEM_malloc_arrayN(
+        singularity_edges_len, sizeof(*singularity_edges), "singularity_edges in solidify");
+    for (NewEdgeRef ***new_edges = orig_edge_data_arr; i < numEdges; i++, new_edges++) {
+      if (*new_edges && (do_shell || edge_adj_faces_len[i] == 1) && (**new_edges)->old_edge == i) {
+        for (NewEdgeRef **l = *new_edges; *l; l++) {
+          if ((*l)->link_edge_groups[0]->is_singularity &&
+              (*l)->link_edge_groups[1]->is_singularity) {
+            const uint v1 = (*l)->link_edge_groups[0]->new_vert;
+            const uint v2 = (*l)->link_edge_groups[1]->new_vert;
+            bool exists_already = false;
+            uint j = 0;
+            for (uint(*p)[2] = singularity_edges; j < totsingularity; p++, j++) {
+              if (((*p)[0] == v1 && (*p)[1] == v2) || ((*p)[0] == v2 && (*p)[1] == v1)) {
+                exists_already = true;
+                break;
+              }
+            }
+            if (!exists_already) {
+              has_singularities = true;
+              if (singularity_edges_len <= totsingularity) {
+                singularity_edges_len = totsingularity + 1;
+                singularity_edges = MEM_reallocN_id(singularity_edges,
+                                                    singularity_edges_len *
+                                                        sizeof(*singularity_edges),
+                                                    "singularity_edges in solidify");
+              }
+              singularity_edges[totsingularity][0] = v1;
+              singularity_edges[totsingularity][1] = v2;
+              totsingularity++;
+              if (edge_adj_faces_len[i] == 1 && do_rim) {
+                numNewLoops -= 2;
+                numNewPolys--;
+              }
+            }
+            else {
+              numNewEdges--;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  /* Create Mesh *result with proper capacity. */
+  result = BKE_mesh_new_nomain_from_template(
+      mesh, (int)(numNewVerts), (int)(numNewEdges), 0, (int)(numNewLoops), (int)(numNewPolys));
+
+  mpoly = result->mpoly;
+  mloop = result->mloop;
+  medge = result->medge;
+  mvert = result->mvert;
+
+  int *origindex_edge = CustomData_get_layer(&result->edata, CD_ORIGINDEX);
+  int *origindex_poly = CustomData_get_layer(&result->pdata, CD_ORIGINDEX);
+
+  /* Make_new_verts. */
+  {
+    gs_ptr = orig_vert_groups_arr;
+    for (uint i = 0; i < numVerts; i++, gs_ptr++) {
+      EdgeGroup *gs = *gs_ptr;
+      if (gs) {
+        EdgeGroup *g = gs;
+        for (uint j = 0; g->valid; j++, g++) {
+          if (g->new_vert != MOD_SOLIDIFY_EMPTY_TAG) {
+            CustomData_copy_data(&mesh->vdata, &result->vdata, (int)i, (int)g->new_vert, 1);
+            copy_v3_v3(mvert[g->new_vert].co, g->co);
+            mvert[g->new_vert].flag = orig_mvert[i].flag;
+          }
+        }
+      }
+    }
+  }
+
+  result->runtime.cd_dirty_vert |= CD_MASK_NORMAL;
+
+  /* Make edges. */
+  {
+    uint i = 0;
+    edge_index += totsingularity;
+    for (NewEdgeRef ***new_edges = orig_edge_data_arr; i < numEdges; i++, new_edges++) {
+      if (*new_edges && (do_shell || edge_adj_faces_len[i] == 1) && (**new_edges)->old_edge == i) {
+        for (NewEdgeRef **l = *new_edges; *l; l++) {
+          if ((*l)->new_edge != MOD_SOLIDIFY_EMPTY_TAG) {
+            const uint v1 = (*l)->link_edge_groups[0]->new_vert;
+            const uint v2 = (*l)->link_edge_groups[1]->new_vert;
+            uint insert = edge_index;
+            if (has_singularities && ((*l)->link_edge_groups[0]->is_singularity &&
+                                      (*l)->link_edge_groups[1]->is_singularity)) {
+              uint j = 0;
+              for (uint(*p)[2] = singularity_edges; j < totsingularity; p++, j++) {
+                if (((*p)[0] == v1 && (*p)[1] == v2) || ((*p)[0] == v2 && (*p)[1] == v1)) {
+                  insert = j;
+                  break;
+                }
+              }
+              BLI_assert(insert == j);
+            }
+            else {
+              edge_index++;
+            }
+            CustomData_copy_data(&mesh->edata, &result->edata, (int)i, (int)insert, 1);
+            BLI_assert(v1 != MOD_SOLIDIFY_EMPTY_TAG);
+            BLI_assert(v2 != MOD_SOLIDIFY_EMPTY_TAG);
+            medge[insert].v1 = v1;
+            medge[insert].v2 = v2;
+            medge[insert].flag = orig_medge[(*l)->old_edge].flag | ME_EDGEDRAW | ME_EDGERENDER;
+            medge[insert].crease = orig_medge[(*l)->old_edge].crease;
+            medge[insert].bweight = orig_medge[(*l)->old_edge].bweight;
+            (*l)->new_edge = insert;
+          }
+        }
+      }
+    }
+  }
+  if (singularity_edges) {
+    MEM_freeN(singularity_edges);
+  }
+
+  /* DEBUG CODE FOR BUGFIXING (can not be removed because every bugfix needs this badly!). */
+#if 0
+  {
+    gs_ptr = orig_vert_groups_arr;
+    for (uint i = 0; i < numVerts; i++, gs_ptr++) {
+      EdgeGroup *gs = *gs_ptr;
+      if (gs) {
+        for (EdgeGroup *g = gs; g->valid; g++) {
+          NewEdgeRef **e = g->edges;
+          for (uint j = 0; j < g->edges_len; j++, e++) {
+            printf("%u/%d, ", (*e)->old_edge, (int)(*e)->new_edge);
+          }
+          printf("(tg:%u)(s:%u,c:%d)\n", g->topo_group, g->split, g->is_orig_closed);
+        }
+        printf("\n");
+      }
+    }
+  }
+#endif
+
+  /* Make boundary edges/faces. */
+  {
+    gs_ptr = orig_vert_groups_arr;
+    for (uint i = 0; i < numVerts; i++, gs_ptr++) {
+      EdgeGroup *gs = *gs_ptr;
+      if (gs) {
+        EdgeGroup *g = gs;
+        EdgeGroup *g2 = gs;
+        EdgeGroup *last_g = NULL;
+        EdgeGroup *first_g = NULL;
+        /* Data calculation cache. */
+        char max_crease;
+        char last_max_crease = 0;
+        char first_max_crease = 0;
+        char max_bweight;
+        char last_max_bweight = 0;
+        char first_max_bweight = 0;
+        short flag;
+        short last_flag = 0;
+        short first_flag = 0;
+        for (uint j = 0; g->valid; g++) {
+          if ((do_rim && !g->is_orig_closed) || (do_shell && g->split)) {
+            max_crease = 0;
+            max_bweight = 0;
+            flag = 0;
+            for (uint k = 1; k < g->edges_len - 1; k++) {
+              ed = orig_medge + g->edges[k]->old_edge;
+              if (ed->crease > max_crease) {
+                max_crease = ed->crease;
+              }
+              if (ed->bweight > max_bweight) {
+                max_bweight = ed->bweight;
+              }
+              flag |= ed->flag;
+            }
+            if (!first_g) {
+              first_g = g;
+              first_max_crease = max_crease;
+              first_max_bweight = max_bweight;
+              first_flag = flag;
+            }
+            else {
+              last_g->open_face_edge = edge_index;
+              CustomData_copy_data(&mesh->edata,
+                                   &result->edata,
+                                   (int)last_g->edges[0]->old_edge,
+                                   (int)edge_index,
+                                   1);
+              if (origindex_edge) {
+                origindex_edge[edge_index] = ORIGINDEX_NONE;
+              }
+              medge[edge_index].v1 = last_g->new_vert;
+              medge[edge_index].v2 = g->new_vert;
+              medge[edge_index].flag = ME_EDGEDRAW | ME_EDGERENDER |
+                                       ((last_flag | flag) & (ME_SEAM | ME_SHARP));
+              medge[edge_index].crease = MAX2(last_max_crease, max_crease);
+              medge[edge_index++].bweight = MAX2(last_max_bweight, max_bweight);
+            }
+            last_g = g;
+            last_max_crease = max_crease;
+            last_max_bweight = max_bweight;
+            last_flag = flag;
+            j++;
+          }
+          if (!(g + 1)->valid || g->topo_group != (g + 1)->topo_group) {
+            if (j == 2) {
+              last_g->open_face_edge = edge_index - 1;
+            }
+            if (j > 2) {
+              CustomData_copy_data(&mesh->edata,
+                                   &result->edata,
+                                   (int)last_g->edges[0]->old_edge,
+                                   (int)edge_index,
+                                   1);
+              if (origindex_edge) {
+                origindex_edge[edge_index] = ORIGINDEX_NONE;
+              }
+              last_g->open_face_edge = edge_index;
+              medge[edge_index].v1 = last_g->new_vert;
+              medge[edge_index].v2 = first_g->new_vert;
+              medge[edge_index].flag = ME_EDGEDRAW | ME_EDGERENDER |
+                                       ((last_flag | first_flag) & (ME_SEAM | ME_SHARP));
+              medge[edge_index].crease = MAX2(last_max_crease, first_max_crease);
+              medge[edge_index++].bweight = MAX2(last_max_bweight, first_max_bweight);
+
+              /* Loop data. */
+              int *loops = MEM_malloc_arrayN(j, sizeof(*loops), "loops in solidify");
+              /* The #mat_nr is from consensus. */
+              short most_mat_nr = 0;
+              uint most_mat_nr_face = 0;
+              uint most_mat_nr_count = 0;
+              for (short l = 0; l < mat_nrs; l++) {
+                uint count = 0;
+                uint face = 0;
+                uint k = 0;
+                for (EdgeGroup *g3 = g2; g3->valid && k < j; g3++) {
+                  if ((do_rim && !g3->is_orig_closed) || (do_shell && g3->split)) {
+                    /* Check both far ends in terms of faces of an edge group. */
+                    if (g3->edges[0]->faces[0]->face->mat_nr == l) {
+                      face = g3->edges[0]->faces[0]->index;
+                      count++;
+                    }
+                    NewEdgeRef *le = g3->edges[g3->edges_len - 1];
+                    if (le->faces[1] && le->faces[1]->face->mat_nr == l) {
+                      face = le->faces[1]->index;
+                      count++;
+                    }
+                    else if (!le->faces[1] && le->faces[0]->face->mat_nr == l) {
+                      face = le->faces[0]->index;
+                      count++;
+                    }
+                    k++;
+                  }
+                }
+                if (count > most_mat_nr_count) {
+                  most_mat_nr = l;
+                  most_mat_nr_face = face;
+                  most_mat_nr_count = count;
+                }
+              }
+              CustomData_copy_data(
+                  &mesh->pdata, &result->pdata, (int)most_mat_nr_face, (int)poly_index, 1);
+              if (origindex_poly) {
+                origindex_poly[poly_index] = ORIGINDEX_NONE;
+              }
+              mpoly[poly_index].loopstart = (int)loop_index;
+              mpoly[poly_index].totloop = (int)j;
+              mpoly[poly_index].mat_nr = most_mat_nr +
+                                         (g->is_orig_closed || !do_rim ? 0 : mat_ofs_rim);
+              CLAMP(mpoly[poly_index].mat_nr, 0, mat_nr_max);
+              mpoly[poly_index].flag = orig_mpoly[most_mat_nr_face].flag;
+              poly_index++;
+
+              for (uint k = 0; g2->valid && k < j; g2++) {
+                if ((do_rim && !g2->is_orig_closed) || (do_shell && g2->split)) {
+                  MPoly *face = g2->edges[0]->faces[0]->face;
+                  ml = orig_mloop + face->loopstart;
+                  for (int l = 0; l < face->totloop; l++, ml++) {
+                    if (vm[ml->v] == i) {
+                      loops[k] = face->loopstart + l;
+                      break;
+                    }
+                  }
+                  k++;
+                }
+              }
+
+              if (!do_flip) {
+                for (uint k = 0; k < j; k++) {
+                  CustomData_copy_data(&mesh->ldata, &result->ldata, loops[k], (int)loop_index, 1);
+                  mloop[loop_index].v = medge[edge_index - j + k].v1;
+                  mloop[loop_index++].e = edge_index - j + k;
+                }
+              }
+              else {
+                for (uint k = 1; k <= j; k++) {
+                  CustomData_copy_data(
+                      &mesh->ldata, &result->ldata, loops[j - k], (int)loop_index, 1);
+                  mloop[loop_index].v = medge[edge_index - k].v2;
+                  mloop[loop_index++].e = edge_index - k;
+                }
+              }
+              MEM_freeN(loops);
+            }
+            /* Reset everything for the next poly. */
+            j = 0;
+            last_g = NULL;
+            first_g = NULL;
+            last_max_crease = 0;
+            first_max_crease = 0;
+            last_max_bweight = 0;
+            first_max_bweight = 0;
+            last_flag = 0;
+            first_flag = 0;
+          }
+        }
+      }
+    }
+  }
+
+  /* Make boundary faces. */
+  if (do_rim) {
+    for (uint i = 0; i < numEdges; i++) {
+      if (edge_adj_faces_len[i] == 1 && orig_edge_data_arr[i] &&
+          (*orig_edge_data_arr[i])->old_edge == i) {
+        NewEdgeRef **new_edges = orig_edge_data_arr[i];
+
+        NewEdgeRef *edge1 = new_edges[0];
+        NewEdgeRef *edge2 = new_edges[1];
+        const bool v1_singularity = edge1->link_edge_groups[0]->is_singularity;
+        const bool v2_singularity = edge1->link_edge_groups[1]->is_singularity;
+        if (v1_singularity && v2_singularity) {
+          continue;
+        }
+
+        MPoly *face = (*new_edges)->faces[0]->face;
+        CustomData_copy_data(
+            &mesh->pdata, &result->pdata, (int)(*new_edges)->faces[0]->index, (int)poly_index, 1);
+        mpoly[poly_index].loopstart = (int)loop_index;
+        mpoly[poly_index].totloop = 4 - (int)(v1_singularity || v2_singularity);
+        mpoly[poly_index].mat_nr = face->mat_nr + mat_ofs_rim;
+        CLAMP(mpoly[poly_index].mat_nr, 0, mat_nr_max);
+        mpoly[poly_index].flag = face->flag;
+        poly_index++;
+
+        int loop1 = -1;
+        int loop2 = -1;
+        ml = orig_mloop + face->loopstart;
+        const uint old_v1 = vm[orig_medge[edge1->old_edge].v1];
+        const uint old_v2 = vm[orig_medge[edge1->old_edge].v2];
+        for (uint j = 0; j < face->totloop; j++, ml++) {
+          if (vm[ml->v] == old_v1) {
+            loop1 = face->loopstart + (int)j;
+          }
+          else if (vm[ml->v] == old_v2) {
+            loop2 = face->loopstart + (int)j;
+          }
+        }
+        BLI_assert(loop1 != -1 && loop2 != -1);
+        MEdge *open_face_edge;
+        uint open_face_edge_index;
+        if (!do_flip) {
+          CustomData_copy_data(&mesh->ldata, &result->ldata, loop1, (int)loop_index, 1);
+          mloop[loop_index].v = medge[edge1->new_edge].v1;
+          mloop[loop_index++].e = edge1->new_edge;
+
+          if (!v2_singularity) {
+            open_face_edge_index = edge1->link_edge_groups[1]->open_face_edge;
+            CustomData_copy_data(&mesh->ldata, &result->ldata, loop2, (int)loop_index, 1);
+            mloop[loop_index].v = medge[edge1->new_edge].v2;
+            open_face_edge = medge + open_face_edge_index;
+            if (ELEM(medge[edge2->new_edge].v2, open_face_edge->v1, open_face_edge->v2)) {
+              mloop[loop_index++].e = open_face_edge_index;
+            }
+            else {
+              mloop[loop_index++].e = edge2->link_edge_groups[1]->open_face_edge;
+            }
+          }
+
+          CustomData_copy_data(&mesh->ldata, &result->ldata, loop2, (int)loop_index, 1);
+          mloop[loop_index].v = medge[edge2->new_edge].v2;
+          mloop[loop_index++].e = edge2->new_edge;
+
+          if (!v1_singularity) {
+            open_face_edge_index = edge2->link_edge_groups[0]->open_face_edge;
+            CustomData_copy_data(&mesh->ldata, &result->ldata, loop1, (int)loop_index, 1);
+            mloop[loop_index].v = medge[edge2->new_edge].v1;
+            open_face_edge = medge + open_face_edge_index;
+            if (ELEM(medge[edge1->new_edge].v1, open_face_edge->v1, open_face_edge->v2)) {
+              mloop[loop_index++].e = open_face_edge_index;
+            }
+            else {
+              mloop[loop_index++].e = edge1->link_edge_groups[0]->open_face_edge;
+            }
+          }
+        }
+        else {
+          if (!v1_singularity) {
+            open_face_edge_index = edge1->link_edge_groups[0]->open_face_edge;
+            CustomData_copy_data(&mesh->ldata, &result->ldata, loop1, (int)loop_index, 1);
+            mloop[loop_index].v = medge[edge1->new_edge].v1;
+            open_face_edge = medge + open_face_edge_index;
+            if (ELEM(medge[edge2->new_edge].v1, open_face_edge->v1, open_face_edge->v2)) {
+              mloop[loop_index++].e = open_face_edge_index;
+            }
+            else {
+              mloop[loop_index++].e = edge2->link_edge_groups[0]->open_face_edge;
+            }
+          }
+
+          CustomData_copy_data(&mesh->ldata, &result->ldata, loop1, (int)loop_index, 1);
+          mloop[loop_index].v = medge[edge2->new_edge].v1;
+          mloop[loop_index++].e = edge2->new_edge;
+
+          if (!v2_singularity) {
+            open_face_edge_index = edge2->link_edge_groups[1]->open_face_edge;
+            CustomData_copy_data(&mesh->ldata, &result->ldata, loop2, (int)loop_index, 1);
+            mloop[loop_index].v = medge[edge2->new_edge].v2;
+            open_face_edge = medge + open_face_edge_index;
+            if (ELEM(medge[edge1->new_edge].v2, open_face_edge->v1, open_face_edge->v2)) {
+              mloop[loop_index++].e = open_face_edge_index;
+            }
+            else {
+              mloop[loop_index++].e = edge1->link_edge_groups[1]->open_face_edge;
+            }
+          }
+
+          CustomData_copy_data(&mesh->ldata, &result->ldata, loop2, (int)loop_index, 1);
+          mloop[loop_index].v = medge[edge1->new_edge].v2;
+          mloop[loop_index++].e = edge1->new_edge;
+        }
+      }
+    }
+  }
+
+  /* Make faces. */
+  if (do_shell) {
+    NewFaceRef *fr = face_sides_arr;
+    uint *face_loops = MEM_malloc_arrayN(
+        largest_ngon * 2, sizeof(*face_loops), "face_loops in solidify");
+    uint *face_verts = MEM_malloc_arrayN(
+        largest_ngon * 2, sizeof(*face_verts), "face_verts in solidify");
+    uint *face_edges = MEM_malloc_arrayN(
+        largest_ngon * 2, sizeof(*face_edges), "face_edges in solidify");
+    for (uint i = 0; i < numPolys * 2; i++, fr++) {
+      const uint loopstart = (uint)fr->face->loopstart;
+      uint totloop = (uint)fr->face->totloop;
+      uint valid_edges = 0;
+      uint k = 0;
+      while (totloop > 0 && (!fr->link_edges[totloop - 1] ||
+                             fr->link_edges[totloop - 1]->new_edge == MOD_SOLIDIFY_EMPTY_TAG)) {
+        totloop--;
+      }
+      if (totloop > 0) {
+        NewEdgeRef *prior_edge = fr->link_edges[totloop - 1];
+        uint prior_flip = (uint)(vm[orig_medge[prior_edge->old_edge].v1] ==
+                                 vm[orig_mloop[loopstart + (totloop - 1)].v]);
+        for (uint j = 0; j < totloop; j++) {
+          NewEdgeRef *new_edge = fr->link_edges[j];
+          if (new_edge && new_edge->new_edge != MOD_SOLIDIFY_EMPTY_TAG) {
+            valid_edges++;
+            const uint flip = (uint)(vm[orig_medge[new_edge->old_edge].v2] ==
+                                     vm[orig_mloop[loopstart + j].v]);
+            BLI_assert(flip ||
+                       vm[orig_medge[new_edge->old_edge].v1] == vm[orig_mloop[loopstart + j].v]);
+            /* The vert thats in the current loop. */
+            const uint new_v1 = new_edge->link_edge_groups[flip]->new_vert;
+            /* The vert thats in the next loop. */
+            const uint new_v2 = new_edge->link_edge_groups[1 - flip]->new_vert;
+            if (k == 0 || face_verts[k - 1] != new_v1) {
+              face_loops[k] = loopstart + j;
+              if (fr->reversed) {
+                face_edges[k] = prior_edge->link_edge_groups[prior_flip]->open_face_edge;
+              }
+              else {
+                face_edges[k] = new_edge->link_edge_groups[flip]->open_face_edge;
+              }
+              BLI_assert(k == 0 || medge[face_edges[k]].v2 == face_verts[k - 1] ||
+                         medge[face_edges[k]].v1 == face_verts[k - 1]);
+              BLI_assert(face_edges[k] == MOD_SOLIDIFY_EMPTY_TAG ||
+                         medge[face_edges[k]].v2 == new_v1 || medge[face_edges[k]].v1 == new_v1);
+              face_verts[k++] = new_v1;
+            }
+            prior_edge = new_edge;
+            prior_flip = 1 - flip;
+            if (j < totloop - 1 || face_verts[0] != new_v2) {
+              face_loops[k] = loopstart + (j + 1) % totloop;
+              face_edges[k] = new_edge->new_edge;
+              face_verts[k++] = new_v2;
+            }
+            else {
+              face_edges[0] = new_edge->new_edge;
+            }
+          }
+        }
+        if (k > 2 && valid_edges > 2) {
+          CustomData_copy_data(&mesh->pdata, &result->pdata, (int)(i / 2), (int)poly_index, 1);
+          mpoly[poly_index].loopstart = (int)loop_index;
+          mpoly[poly_index].totloop = (int)k;
+          mpoly[poly_index].mat_nr = fr->face->mat_nr + mat_ofs;
+          CLAMP(mpoly[poly_index].mat_nr, 0, mat_nr_max);
+          mpoly[poly_index].flag = fr->face->flag;
+          if (fr->reversed != do_flip) {
+            for (int l = (int)k - 1; l >= 0; l--) {
+              CustomData_copy_data(
+                  &mesh->ldata, &result->ldata, (int)face_loops[l], (int)loop_index, 1);
+              mloop[loop_index].v = face_verts[l];
+              mloop[loop_index++].e = face_edges[l];
+            }
+          }
+          else {
+            uint l = k - 1;
+            for (uint next_l = 0; next_l < k; next_l++) {
+              CustomData_copy_data(
+                  &mesh->ldata, &result->ldata, (int)face_loops[l], (int)loop_index, 1);
+              mloop[loop_index].v = face_verts[l];
+              mloop[loop_index++].e = face_edges[next_l];
+              l = next_l;
+            }
+          }
+          poly_index++;
+        }
+      }
+    }
+    MEM_freeN(face_loops);
+    MEM_freeN(face_verts);
+    MEM_freeN(face_edges);
+  }
+  if (edge_index != numNewEdges) {
+    modifier_setError(
+        md, "Internal Error: edges array wrong size: %u instead of %u", numNewEdges, edge_index);
+  }
+  if (poly_index != numNewPolys) {
+    modifier_setError(
+        md, "Internal Error: polys array wrong size: %u instead of %u", numNewPolys, poly_index);
+  }
+  if (loop_index != numNewLoops) {
+    modifier_setError(
+        md, "Internal Error: loops array wrong size: %u instead of %u", numNewLoops, loop_index);
+  }
+  BLI_assert(edge_index == numNewEdges);
+  BLI_assert(poly_index == numNewPolys);
+  BLI_assert(loop_index == numNewLoops);
+
+  /* Free remaining memory */
+  {
+    MEM_freeN(vm);
+    MEM_freeN(edge_adj_faces_len);
+    uint i = 0;
+    for (EdgeGroup **p = orig_vert_groups_arr; i < numVerts; i++, p++) {
+      if (*p) {
+        for (EdgeGroup *eg = *p; eg->valid; eg++) {
+          MEM_freeN(eg->edges);
+        }
+        MEM_freeN(*p);
+      }
+    }
+    MEM_freeN(orig_vert_groups_arr);
+    i = numEdges;
+    for (NewEdgeRef ***p = orig_edge_data_arr + (numEdges - 1); i > 0; i--, p--) {
+      if (*p && (**p)->old_edge == i - 1) {
+        for (NewEdgeRef **l = *p; *l; l++) {
+          MEM_freeN(*l);
+        }
+        MEM_freeN(*p);
+      }
+    }
+    MEM_freeN(orig_edge_data_arr);
+    MEM_freeN(orig_edge_lengths);
+    i = 0;
+    for (NewFaceRef *p = face_sides_arr; i < numPolys * 2; i++, p++) {
+      MEM_freeN(p->link_edges);
+    }
+    MEM_freeN(face_sides_arr);
+    MEM_freeN(poly_nors);
+  }
+
+#undef MOD_SOLIDIFY_EMPTY_TAG
+
+  return result;
+}
+
+/** \} */
diff --git a/source/blender/modifiers/intern/MOD_solidify_util.h b/source/blender/modifiers/intern/MOD_solidify_util.h
new file mode 100644
index 00000000000..dba360dc1ce
--- /dev/null
+++ b/source/blender/modifiers/intern/MOD_solidify_util.h
@@ -0,0 +1,34 @@
+/*
+ * 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.
+ */
+
+/** \file
+ * \ingroup modifiers
+ */
+
+#ifndef __MOD_MESHCACHE_UTIL_H__
+#define __MOD_MESHCACHE_UTIL_H__
+
+/* MOD_solidify_extrude.c */
+Mesh *MOD_solidify_extrude_applyModifier(ModifierData *md,
+                                         const ModifierEvalContext *ctx,
+                                         Mesh *mesh);
+
+/* MOD_solidify_nonmanifold.c */
+Mesh *MOD_solidify_nonmanifold_applyModifier(ModifierData *md,
+                                             const ModifierEvalContext *ctx,
+                                             Mesh *mesh);
+
+#endif /* __MOD_MESHCACHE_UTIL_H__ */