profiling

Some of these could be enabled if there is an interest to consider
them warnings, until then they're noisy and don't often hint at errors
in the code.

@@ -1497,10 +1480,12 @@ static void sculpt_expand_update_for_vertex(bContext *C, Object *ob, const PBVHV
           sculpt_expand_mask_update_task(ss, mask_write, expand_cache->nodes[i]);
         }
       });
+      SCULPT_flush_update_step(C, SCULPT_UPDATE_MASK);

@@ -2140,6 +2124,11 @@ static int sculpt_expand_invoke(bContext *C, wmOperator *op, const wmEvent *even
   if (!BKE_base_is_visible(v3d, base)) {
     return OPERATOR_CANCELLED;
   }
+  /* Do nothing when the mesh has 0 vertices. */

@@ -2222,7 +2204,6 @@ static int sculpt_expand_invoke(bContext *C, wmOperator *op, const wmEvent *even
       ss->expand_cache, ob, ss->expand_cache->initial_active_vertex, falloff_type);
 
   sculpt_expand_check_topology_islands(ob, falloff_type);

@@ -8,6 +8,7 @@
 
 #include <cmath>
 #include <cstdlib>
+#include <thread>

@@ -77,6 +78,7 @@ static bool sculpt_geodesic_mesh_test_dist_add(Span<float3> vert_positions,
   return false;
 }
 
+#define UNASSIGNED -1

@@ -146,0 +137,4 @@
+      threading::parallel_for(IndexRange(0, totvert), 4096, [&](IndexRange range) {
+        for (const int i : range) {
+          if (len_squared_v3v3(v_co, vert_positions[i]) <= limit_radius_sq) {
+            affected_vert[i].set();

@@ -147,0 +144,4 @@
+    }
+  }
+
+  std::vector<int> queue, queue_next;

@@ -147,0 +148,4 @@
+  queue.reserve(totedge);
+
+  /* Add edges adjacent to an initial vertex to the queue.
+   * Since initial vertex are few only, iterating over its neighbour edges

@@ -147,0 +149,4 @@
+
+  /* Add edges adjacent to an initial vertex to the queue.
+   * Since initial vertex are few only, iterating over its neighbour edges
+   * instead of over all edges scales better as mesh edge count increases */

@@ -147,0 +153,4 @@
+  GSetIterator gs_iter;
+  GSET_ITER (gs_iter, initial_verts) {
+    const int seed_vert = POINTER_AS_INT(BLI_gsetIterator_getKey(&gs_iter));
+    for (const int e : ss->vert_to_edge_map[seed_vert]) {

@@ -180,3 +191,3 @@
         }
         for (const int v_other : corner_verts.slice(faces[face])) {
-          if (ELEM(v_other, v1, v2)) {
+          if (ELEM(v_other, v1, v2) || !sculpt_geodesic_mesh_test_dist_add(

@@ -202,2 +199,2 @@
-                  BLI_LINKSTACK_PUSH(queue_next, POINTER_FROM_INT(e_other));
-                }
+          for (const int e_other : ss->vert_to_edge_map[v_other]) {
+            const int ev_other = (edges[e_other][0] == v_other) ? edges[e_other][1] :

@@ -204,0 +213,4 @@
+             * collission in hash tables can also be used to avoid
+             * parallel data storage collisions */
+            size_t idx = e_other % new_size;
+            while (queue_next[idx] != UNASSIGNED) {

@@ -303,3 +321,3 @@
 }
 
-}  // namespace blender::ed::sculpt_paint::geodesic
+} /* namespace blender::ed::sculpt_paint::geodesic */

@@ -1 +1 @@
-Subproject commit 612e310bf7ee505a79dad7ea925916ec9a812995
+Subproject commit ce850ef056c57d0e4127eda48a514b2c4a651435

@@ -1 +1 @@
-Subproject commit 612e310bf7ee505a79dad7ea925916ec9a812995
+Subproject commit 42979e6ec67f657c7efce9cb672fe07b3f8fff9d

@@ -610,3 +610,2 @@
       const float *co = SCULPT_vertex_co_get(ss, symm_vertex);
-      for (int i = 0; i < totvert; i++) {
-        PBVHVertRef vertex = BKE_pbvh_index_to_vertex(ss->pbvh, i);
+      threading::parallel_for(IndexRange(0, totvert), 4096, [&](IndexRange range) {

@@ -79,14 +79,13 @@ static bool sculpt_geodesic_mesh_test_dist_add(Span<float3> vert_positions,
 
 static float *geodesic_mesh_create(Object *ob, GSet *initial_verts, const float limit_radius)
 {
+  const int UNASSIGNED = -1;

@@ -120,0 +111,4 @@
+        dists[i] = 0.0f;
+      }
+      else {
+        dists[i] = FLT_MAX;

@@ -147,0 +143,4 @@
+    }
+  }
+
+  Vector<int> queue, queue_next;

@@ -161,1 +167,3 @@
-  }
+  /* queue and queue_next should not differ by a huge amount of size
+   * since they are an advancing front in the mesh hence no need
+   * to allocate much more memmory and keep the iteration range smaller */

@@ -167,1 +173,3 @@
-      int v2 = edges[e][1];
+  BitVector<> edge_tag(totedge);
+  while (!queue.is_empty()) {
+    threading::parallel_for_each(IndexRange(0, queue.size()), [&](const int val) {

@@ -168,0 +174,4 @@
+  while (!queue.is_empty()) {
+    threading::parallel_for_each(IndexRange(0, queue.size()), [&](const int val) {
+      const int edge_index = queue[val];
+      int v1 = edges[edge_index][0];

@@ -204,0 +216,4 @@
+            /* Open addressing with linear probing that minimizes
+             * collission in hash tables can also be used to avoid
+             * parallel data storage collisions */
+            size_t idx = edge_other % new_size;

@@ -80,0 +80,4 @@
+/* Open addressing with linear probing that minimizes
+ * collission in hash tables can also be used to avoid
+ * parallel data storage collisions */
+static void lock_free_emplace(const int edge_other, const int non_checked, Vector<int> &queue_next)

@@ -147,0 +160,4 @@
+  }
+
+  Vector<int> queue;
+  Vector<int> queue_next;

@@ -209,0 +245,4 @@
+    });
+
+    queue.clear();
+    for (int i = 0; i < new_size; ++i) {

@@ -80,0 +92,4 @@
+      idx = 0;
+    }
+  }
+  queue_next[idx] = edge_other;

@@ -172,0 +204,4 @@
+          if (dists[edge_vert_a] > dists[edge_vert_b]) {
+            std::swap(edge_vert_a, edge_vert_b);
+          }
+          sculpt_geodesic_mesh_test_dist_add(vert_positions,

@@ -2141,6 +2143,11 @@ static int sculpt_expand_invoke(bContext *C, wmOperator *op, const wmEvent *even
     return OPERATOR_CANCELLED;
   }
 
+  const int totvert = SCULPT_vertex_count_get(ss);

@@ -50,3 +50,3 @@
   }
 
-  BLI_assert(dists[v1] != FLT_MAX);
+  BLI_assert(dists[v1] != std::numeric_limits<float>::max());

@@ -77,16 +77,33 @@ static bool sculpt_geodesic_mesh_test_dist_add(Span<float3> vert_positions,
   return false;
 }
 
+/* Open addressing with linear probing that minimizes

@@ -161,1 +186,3 @@
-  }
+  /* queue and queue_next should not differ by a huge amount of size
+   * since they are an advancing front in the mesh hence no need
+   * to allocate much more memmory and keep the iteration range smaller */

@@ -52,2 +65,2 @@
-  BLI_assert(dists[v1] != FLT_MAX);
-  if (dists[v0] <= dists[v1]) {
+  BLI_assert(dists[v1].load() != std::numeric_limits<float>::max());
+  if (dists[v0].load() <= dists[v1].load()) {

@@ -74,0 +80,4 @@
+     * the actual order of update is irrelevant here */
+    std::atomic<float> &atomic_dist_v0 = reinterpret_cast<std::atomic<float> &>(dists[v0]);
+    if (atomic_dist_v0.load() > dist0) {
+      atomic_dist_v0.store(dist0);