release/datafiles/userdef/userdef_default_theme.c

@@ -869,6 +869,7 @@ const bTheme U_theme_default = {
     .nodeclass_layout = RGBA(0x6c696fff),
     .nodeclass_geometry = RGBA(0x00d6a3ff),
     .nodeclass_attribute = RGBA(0x001566ff),
+    .node_zone_simulation = RGBA(0x66416233),
     .movie = RGBA(0x0f0f0fcc),
     .gp_vertex_size = 3,
     .gp_vertex = RGBA(0x97979700),

scripts/startup/bl_operators/geometry_nodes.py

@@ -5,6 +5,10 @@ from bpy.types import Operator
 
 from bpy.app.translations import pgettext_data as data_
 
+from bpy.props import (
+    EnumProperty,
+)
+
 
 def build_default_empty_geometry_node_group(name):
     group = bpy.data.node_groups.new(name, 'GeometryNodeTree')
@@ -242,8 +246,102 @@ class NewGeometryNodeTreeAssign(Operator):
         return {'FINISHED'}
 
 
+class SimulationZoneOperator:
+    input_node_type = 'GeometryNodeSimulationInput'
+    output_node_type = 'GeometryNodeSimulationOutput'
+
+    @classmethod
+    def get_output_node(cls, context):
+        node = context.active_node
+        if node.bl_idname == cls.input_node_type:
+            return node.paired_output
+        if node.bl_idname == cls.output_node_type:
+            return node
+
+    @classmethod
+    def poll(cls, context):
+        space = context.space_data
+        # Needs active node editor and a tree.
+        if not space or space.type != 'NODE_EDITOR' or not space.edit_tree or space.edit_tree.library:
+            return False
+        node = context.active_node
+        if node is None or node.bl_idname not in [cls.input_node_type, cls.output_node_type]:
+            return False
+        if cls.get_output_node(context) is None:
+            return False
+        return True
+
+
+class SimulationZoneItemAddOperator(SimulationZoneOperator, Operator):
+    """Add a state item to the simulation zone"""
+    bl_idname = "node.simulation_zone_item_add"
+    bl_label = "Add State Item"
+    bl_options = {'REGISTER', 'UNDO'}
+
+    default_socket_type = 'GEOMETRY'
+
+    def execute(self, context):
+        node = self.get_output_node(context)
+        state_items = node.state_items
+
+        # Remember index to move the item.
+        dst_index = min(node.active_index + 1, len(state_items))
+        # Empty name so it is based on the type only.
+        state_items.new(self.default_socket_type, "")
+        state_items.move(len(state_items) - 1, dst_index)
+        node.active_index = dst_index
+
+        return {'FINISHED'}
+
+
+class SimulationZoneItemRemoveOperator(SimulationZoneOperator, Operator):
+    """Remove a state item from the simulation zone"""
+    bl_idname = "node.simulation_zone_item_remove"
+    bl_label = "Remove State Item"
+    bl_options = {'REGISTER', 'UNDO'}
+
+    def execute(self, context):
+        node = self.get_output_node(context)
+        state_items = node.state_items
+
+        if node.active_item:
+            state_items.remove(node.active_item)
+            node.active_index = min(node.active_index, len(state_items) - 1)
+
+        return {'FINISHED'}
+
+
+class SimulationZoneItemMoveOperator(SimulationZoneOperator, Operator):
+    """Move a simulation state item up or down in the list"""
+    bl_idname = "node.simulation_zone_item_move"
+    bl_label = "Move State Item"
+    bl_options = {'REGISTER', 'UNDO'}
+
+    direction: EnumProperty(
+        name="Direction",
+        items=[('UP', "Up", ""), ('DOWN', "Down", "")],
+        default='UP',
+    )
+
+    def execute(self, context):
+        node = self.get_output_node(context)
+        state_items = node.state_items
+
+        if self.direction == 'UP' and node.active_index > 0:
+            state_items.move(node.active_index, node.active_index - 1)
+            node.active_index = node.active_index - 1
+        elif self.direction == 'DOWN' and node.active_index < len(state_items) - 1:
+            state_items.move(node.active_index, node.active_index + 1)
+            node.active_index = node.active_index + 1
+
+        return {'FINISHED'}
+
+
 classes = (
     NewGeometryNodesModifier,
     NewGeometryNodeTreeAssign,
     MoveModifierToNodes,
+    SimulationZoneItemAddOperator,
+    SimulationZoneItemRemoveOperator,
+    SimulationZoneItemMoveOperator,
 )

scripts/startup/bl_operators/node.py

@@ -9,8 +9,12 @@ from bpy.types import (
 from bpy.props import (
     BoolProperty,
     CollectionProperty,
+    FloatVectorProperty,
     StringProperty,
 )
+from mathutils import (
+    Vector,
+)
 
 from bpy.app.translations import pgettext_tip as tip_
 
@@ -27,10 +31,6 @@ class NodeSetting(PropertyGroup):
 # Base class for node "Add" operators.
 class NodeAddOperator:
 
-    type: StringProperty(
-        name="Node Type",
-        description="Node type",
-    )
     use_transform: BoolProperty(
         name="Use Transform",
         description="Start transform operator after inserting the node",
@@ -56,21 +56,18 @@ class NodeAddOperator:
         else:
             space.cursor_location = tree.view_center
 
-    # XXX explicit node_type argument is usually not necessary,
-    # but required to make search operator work:
-    # add_search has to override the 'type' property
-    # since it's hardcoded in bpy_operator_wrap.c ...
-    def create_node(self, context, node_type=None):
+    # Deselect all nodes in the tree.
+    @staticmethod
+    def deselect_nodes(context):
         space = context.space_data
         tree = space.edit_tree
-
-        if node_type is None:
-            node_type = self.type
-
-        # select only the new node
         for n in tree.nodes:
             n.select = False
 
+    def create_node(self, context, node_type):
+        space = context.space_data
+        tree = space.edit_tree
+
         try:
             node = tree.nodes.new(type=node_type)
         except RuntimeError as e:
@@ -109,14 +106,6 @@ class NodeAddOperator:
         return (space and (space.type == 'NODE_EDITOR') and
                 space.edit_tree and not space.edit_tree.library)
 
-    # Default execute simply adds a node
-    def execute(self, context):
-        if self.properties.is_property_set("type"):
-            self.create_node(context)
-            return {'FINISHED'}
-        else:
-            return {'CANCELLED'}
-
     # Default invoke stores the mouse position to place the node correctly
     # and optionally invokes the transform operator
     def invoke(self, context, event):
@@ -129,6 +118,28 @@ class NodeAddOperator:
 
         return result
 
+
+# Simple basic operator for adding a node.
+class NODE_OT_add_node(NodeAddOperator, Operator):
+    """Add a node to the active tree"""
+    bl_idname = "node.add_node"
+    bl_label = "Add Node"
+    bl_options = {'REGISTER', 'UNDO'}
+
+    type: StringProperty(
+        name="Node Type",
+        description="Node type",
+    )
+
+    # Default execute simply adds a node.
+    def execute(self, context):
+        if self.properties.is_property_set("type"):
+            self.deselect_nodes(context)
+            self.create_node(context, self.type)
+            return {'FINISHED'}
+        else:
+            return {'CANCELLED'}
+
     @classmethod
     def description(cls, _context, properties):
         nodetype = properties["type"]
@@ -139,13 +150,46 @@ class NodeAddOperator:
             return ""
 
 
-# Simple basic operator for adding a node
-class NODE_OT_add_node(NodeAddOperator, Operator):
-    """Add a node to the active tree"""
-    bl_idname = "node.add_node"
-    bl_label = "Add Node"
+class NODE_OT_add_simulation_zone(NodeAddOperator, Operator):
+    """Add simulation zone input and output nodes to the active tree"""
+    bl_idname = "node.add_simulation_zone"
+    bl_label = "Add Simulation Zone"
     bl_options = {'REGISTER', 'UNDO'}
 
+    input_node_type = "GeometryNodeSimulationInput"
+    output_node_type = "GeometryNodeSimulationOutput"
+    offset: FloatVectorProperty(
+        name="Offset",
+        description="Offset of nodes from the cursor when added",
+        size=2,
+        default=(150, 0),
+    )
+
+    def execute(self, context):
+        space = context.space_data
+        tree = space.edit_tree
+
+        props = self.properties
+
+        self.deselect_nodes(context)
+        input_node = self.create_node(context, self.input_node_type)
+        output_node = self.create_node(context, self.output_node_type)
+        if input_node is None or output_node is None:
+            return {'CANCELLED'}
+
+        # Simulation input must be paired with the output.
+        input_node.pair_with_output(output_node)
+
+        input_node.location -= Vector(self.offset)
+        output_node.location += Vector(self.offset)
+
+        # Connect geometry sockets by default.
+        from_socket = input_node.outputs.get("Geometry")
+        to_socket = output_node.inputs.get("Geometry")
+        tree.links.new(to_socket, from_socket)
+
+        return {'FINISHED'}
+
 
 class NODE_OT_collapse_hide_unused_toggle(Operator):
     """Toggle collapsed nodes and hide unused sockets"""
@@ -202,6 +246,7 @@ classes = (
     NodeSetting,
 
     NODE_OT_add_node,
+    NODE_OT_add_simulation_zone,
     NODE_OT_collapse_hide_unused_toggle,
     NODE_OT_tree_path_parent,
 )

scripts/startup/bl_ui/__init__.py

@@ -41,6 +41,7 @@ _modules = [
     "properties_physics_common",
     "properties_physics_dynamicpaint",
     "properties_physics_field",
+    "properties_physics_geometry_nodes",
     "properties_physics_rigidbody",
     "properties_physics_rigidbody_constraint",
     "properties_physics_fluid",

scripts/startup/bl_ui/node_add_menu.py

@@ -57,6 +57,18 @@ def draw_root_assets(layout):
     layout.menu_contents("NODE_MT_node_add_root_catalogs")
 
 
+def add_simulation_zone(layout, label):
+    """Add simulation zone to a menu."""
+    target_bl_rna = bpy.types.Node.bl_rna_get_subclass("GeometryNodeSimulationOutput")
+    if target_bl_rna:
+        translation_context = target_bl_rna.translation_context
+    else:
+        translation_context = i18n_contexts.default
+    props = layout.operator("node.add_simulation_zone", text=label, text_ctxt=translation_context)
+    props.use_transform = True
+    return props
+
+
 classes = (
 )
 

scripts/startup/bl_ui/node_add_menu_geometry.py

@@ -464,6 +464,16 @@ class NODE_MT_category_GEO_POINT(Menu):
         node_add_menu.draw_assets_for_catalog(layout, self.bl_label)
 
 
+class NODE_MT_category_simulation(Menu):
+    bl_idname = "NODE_MT_category_simulation"
+    bl_label = "Simulation"
+
+    def draw(self, _context):
+        layout = self.layout
+        node_add_menu.add_simulation_zone(layout, label="Simulation Zone")
+        node_add_menu.draw_assets_for_catalog(layout, self.bl_label)
+
+
 class NODE_MT_category_GEO_TEXT(Menu):
     bl_idname = "NODE_MT_category_GEO_TEXT"
     bl_label = "Text"
@@ -646,6 +656,8 @@ class NODE_MT_geometry_node_add_all(Menu):
         layout.menu("NODE_MT_category_GEO_POINT")
         layout.menu("NODE_MT_category_GEO_VOLUME")
         layout.separator()
+        layout.menu("NODE_MT_category_simulation")
+        layout.separator()
         layout.menu("NODE_MT_geometry_node_GEO_MATERIAL")
         layout.menu("NODE_MT_category_GEO_TEXTURE")
         layout.menu("NODE_MT_category_GEO_UTILITIES")
@@ -685,6 +697,7 @@ classes = (
     NODE_MT_category_PRIMITIVES_MESH,
     NODE_MT_geometry_node_mesh_topology,
     NODE_MT_category_GEO_POINT,
+    NODE_MT_category_simulation,
     NODE_MT_category_GEO_VOLUME,
     NODE_MT_geometry_node_GEO_MATERIAL,
     NODE_MT_category_GEO_TEXTURE,

scripts/startup/bl_ui/properties_physics_geometry_nodes.py

@@ -0,0 +1,47 @@
+# SPDX-License-Identifier: GPL-2.0-or-later
+
+from bpy.types import (
+    Panel,
+)
+from bpy.app.translations import pgettext_iface as iface_
+
+
+class PHYSICS_PT_geometry_nodes(Panel):
+    bl_space_type = 'PROPERTIES'
+    bl_region_type = 'WINDOW'
+    bl_context = "physics"
+    bl_label = "Simulation Nodes"
+    bl_options = {'DEFAULT_CLOSED'}
+
+    @classmethod
+    def geometry_nodes_objects(cls, context):
+        for ob in context.selected_editable_objects:
+            if any([modifier.type == 'NODES' for modifier in ob.modifiers]):
+                yield ob
+
+    @classmethod
+    def poll(cls, context):
+        return any(cls.geometry_nodes_objects(context))
+
+    def draw(self, context):
+        layout = self.layout
+
+        if len(context.selected_editable_objects) > 1:
+            bake_text = iface_("Bake Selected")
+        else:
+            bake_text = iface_("Bake")
+
+        row = layout.row(align=True)
+        row.operator("object.simulation_nodes_cache_bake", text=bake_text).selected = True
+        row.operator("object.simulation_nodes_cache_delete", text="", icon='TRASH').selected = True
+
+
+classes = (
+    PHYSICS_PT_geometry_nodes,
+)
+
+
+if __name__ == "__main__":  # only for live edit.
+    from bpy.utils import register_class
+    for cls in classes:
+        register_class(cls)

scripts/startup/bl_ui/space_node.py

@@ -955,6 +955,85 @@ class NODE_PT_node_tree_interface_outputs(NodeTreeInterfacePanel):
         self.draw_socket_list(context, "OUT", "outputs", "active_output")
 
 
+class NODE_UL_simulation_zone_items(bpy.types.UIList):
+    def draw_item(self, context, layout, _data, item, icon, _active_data, _active_propname, _index):
+        if self.layout_type in {'DEFAULT', 'COMPACT'}:
+            row = layout.row(align=True)
+
+            row.template_node_socket(color=item.color)
+            row.prop(item, "name", text="", emboss=False, icon_value=icon)
+        elif self.layout_type == 'GRID':
+            layout.alignment = 'CENTER'
+            layout.template_node_socket(color=item.color)
+
+
+class NODE_PT_simulation_zone_items(Panel):
+    bl_space_type = 'NODE_EDITOR'
+    bl_region_type = 'UI'
+    bl_category = "Node"
+    bl_label = "Simulation State"
+
+    input_node_type = 'GeometryNodeSimulationInput'
+    output_node_type = 'GeometryNodeSimulationOutput'
+
+    @classmethod
+    def get_output_node(cls, context):
+        node = context.active_node
+        if node.bl_idname == cls.input_node_type:
+            return node.paired_output
+        if node.bl_idname == cls.output_node_type:
+            return node
+
+    @classmethod
+    def poll(cls, context):
+        snode = context.space_data
+        if snode is None:
+            return False
+        node = context.active_node
+        if node is None or node.bl_idname not in [cls.input_node_type, cls.output_node_type]:
+            return False
+        if cls.get_output_node(context) is None:
+            return False
+        return True
+
+    def draw(self, context):
+        layout = self.layout
+
+        output_node = self.get_output_node(context)
+
+        split = layout.row()
+
+        split.template_list(
+            "NODE_UL_simulation_zone_items",
+            "",
+            output_node,
+            "state_items",
+            output_node,
+            "active_index")
+
+        ops_col = split.column()
+
+        add_remove_col = ops_col.column(align=True)
+        add_remove_col.operator("node.simulation_zone_item_add", icon='ADD', text="")
+        add_remove_col.operator("node.simulation_zone_item_remove", icon='REMOVE', text="")
+
+        ops_col.separator()
+
+        up_down_col = ops_col.column(align=True)
+        props = up_down_col.operator("node.simulation_zone_item_move", icon='TRIA_UP', text="")
+        props.direction = 'UP'
+        props = up_down_col.operator("node.simulation_zone_item_move", icon='TRIA_DOWN', text="")
+        props.direction = 'DOWN'
+
+        active_item = output_node.active_item
+        if active_item is not None:
+            layout.use_property_split = True
+            layout.use_property_decorate = False
+            layout.prop(active_item, "socket_type")
+            if active_item.socket_type in {'VECTOR', 'INT', 'BOOLEAN', 'FLOAT', 'RGBA'}:
+                layout.prop(active_item, "attribute_domain")
+
+
 # Grease Pencil properties
 class NODE_PT_annotation(AnnotationDataPanel, Panel):
     bl_space_type = 'NODE_EDITOR'
@@ -1019,6 +1098,8 @@ classes = (
     NODE_UL_interface_sockets,
     NODE_PT_node_tree_interface_inputs,
     NODE_PT_node_tree_interface_outputs,
+    NODE_UL_simulation_zone_items,
+    NODE_PT_simulation_zone_items,
 
     node_panel(EEVEE_MATERIAL_PT_settings),
     node_panel(MATERIAL_PT_viewport),

source/blender/blenkernel/BKE_geometry_set.hh

@@ -191,6 +191,11 @@ struct GeometrySet {
    * access to their data, which might be freed later if this geometry set outlasts the data.
    */
   void ensure_owns_direct_data();
+  /**
+   * Same as #ensure_owns_direct_data but also turns object/collection instances into geometry
+   * instances so that they can be owned.
+   */
+  void ensure_owns_all_data();
 
   using AttributeForeachCallback =
       blender::FunctionRef<void(const blender::bke::AttributeIDRef &attribute_id,

source/blender/blenkernel/BKE_node.h

@@ -1584,6 +1584,9 @@ void BKE_nodetree_remove_layer_n(struct bNodeTree *ntree, struct Scene *scene, i
 #define GEO_NODE_MEAN_FILTER_SDF_VOLUME 1197
 #define GEO_NODE_OFFSET_SDF_VOLUME 1198
 #define GEO_NODE_INDEX_OF_NEAREST 1199
+/* Function nodes use the range starting at 1200. */
+#define GEO_NODE_SIMULATION_INPUT 2100
+#define GEO_NODE_SIMULATION_OUTPUT 2101
 
 /** \} */
 

source/blender/blenkernel/BKE_node_runtime.hh

@@ -31,6 +31,9 @@ struct RelationsInNode;
 }
 namespace aal = anonymous_attribute_lifetime;
 }  // namespace blender::nodes
+namespace blender::bke::node_tree_zones {
+class TreeZones;
+}
 
 namespace blender {
 
@@ -64,6 +67,8 @@ struct NodeIDEquality {
 
 namespace blender::bke {
 
+using NodeIDVectorSet = VectorSet<bNode *, DefaultProbingStrategy, NodeIDHash, NodeIDEquality>;
+
 class bNodeTreeRuntime : NonCopyable, NonMovable {
  public:
   /**
@@ -89,7 +94,7 @@ class bNodeTreeRuntime : NonCopyable, NonMovable {
    * allow simpler and more cache friendly iteration. Supports lookup by integer or by node.
    * Unlike other caches, this is maintained eagerly while changing the tree.
    */
-  VectorSet<bNode *, DefaultProbingStrategy, NodeIDHash, NodeIDEquality> nodes_by_id;
+  NodeIDVectorSet nodes_by_id;
 
   /** Execution data.
    *
@@ -136,6 +141,9 @@ class bNodeTreeRuntime : NonCopyable, NonMovable {
    */
   mutable std::atomic<int> allow_use_dirty_topology_cache = 0;
 
+  CacheMutex tree_zones_cache_mutex;
+  std::unique_ptr<node_tree_zones::TreeZones> tree_zones;
+
   /** Only valid when #topology_cache_is_dirty is false. */
   Vector<bNodeLink *> links;
   Vector<bNodeSocket *> sockets;
@@ -361,12 +369,14 @@ inline blender::Span<bNode *> bNodeTree::all_nodes()
 
 inline bNode *bNodeTree::node_by_id(const int32_t identifier)
 {
+  BLI_assert(identifier >= 0);
   bNode *const *node = this->runtime->nodes_by_id.lookup_key_ptr_as(identifier);
   return node ? *node : nullptr;
 }
 
 inline const bNode *bNodeTree::node_by_id(const int32_t identifier) const
 {
+  BLI_assert(identifier >= 0);
   const bNode *const *node = this->runtime->nodes_by_id.lookup_key_ptr_as(identifier);
   return node ? *node : nullptr;
 }

source/blender/blenkernel/BKE_node_tree_zones.hh

@@ -0,0 +1,43 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+/** \file
+ * \ingroup bke
+ */
+
+#include "DNA_node_types.h"
+
+#include "BLI_vector.hh"
+
+namespace blender::bke::node_tree_zones {
+
+struct TreeZone {
+  TreeZones *owner = nullptr;
+  /** Index of the zone in the array of all zones in a node tree. */
+  int index = -1;
+  /** Zero for top level zones, one for a nested zone, and so on. */
+  int depth = -1;
+  /** Input node of the zone. */
+  const bNode *input_node = nullptr;
+  /** Output node of the zone. */
+  const bNode *output_node = nullptr;
+  /** Direct parent of the zone. If this is null, this is a top level zone. */
+  TreeZone *parent_zone = nullptr;
+  /** Direct children zones. Does not contain recursively nested zones. */
+  Vector<TreeZone *> child_zones;
+  /** Direct children nodes. Does not contain recursively nested nodes. */
+  Vector<const bNode *> child_nodes;
+
+  bool contains_node_recursively(const bNode &node) const;
+};
+
+class TreeZones {
+ public:
+  Vector<std::unique_ptr<TreeZone>> zones;
+  Map<int, int> parent_zone_by_node_id;
+};
+
+const TreeZones *get_tree_zones(const bNodeTree &tree);
+
+}  // namespace blender::bke::node_tree_zones

source/blender/blenkernel/BKE_simulation_state.hh

@@ -0,0 +1,168 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "BKE_geometry_set.hh"
+
+#include "BLI_map.hh"
+#include "BLI_sub_frame.hh"
+
+namespace blender::bke::sim {
+
+class BDataSharing;
+class ModifierSimulationCache;
+
+class SimulationStateItem {
+ public:
+  virtual ~SimulationStateItem() = default;
+};
+
+class GeometrySimulationStateItem : public SimulationStateItem {
+ private:
+  GeometrySet geometry_;
+
+ public:
+  GeometrySimulationStateItem(GeometrySet geometry);
+
+  const GeometrySet &geometry() const
+  {
+    return geometry_;
+  }
+
+  GeometrySet &geometry()
+  {
+    return geometry_;
+  }
+};
+
+class AttributeSimulationStateItem : public SimulationStateItem {
+ private:
+  std::string name_;
+
+ public:
+  AttributeSimulationStateItem(std::string name) : name_(std::move(name)) {}
+
+  StringRefNull name() const
+  {
+    return name_;
+  }
+};
+
+class PrimitiveSimulationStateItem : public SimulationStateItem {
+ private:
+  const CPPType &type_;
+  void *value_;
+
+ public:
+  PrimitiveSimulationStateItem(const CPPType &type, const void *value);
+  ~PrimitiveSimulationStateItem();
+
+  const void *value() const
+  {
+    return value_;
+  }
+
+  const CPPType &type() const
+  {
+    return type_;
+  }
+};
+
+class StringSimulationStateItem : public SimulationStateItem {
+ private:
+  std::string value_;
+
+ public:
+  StringSimulationStateItem(std::string value);
+
+  StringRefNull value() const
+  {
+    return value_;
+  }
+};
+
+class SimulationZoneState {
+ public:
+  Map<int, std::unique_ptr<SimulationStateItem>> item_by_identifier;
+};
+
+struct SimulationZoneID {
+  Vector<int> node_ids;
+
+  uint64_t hash() const
+  {
+    return get_default_hash(this->node_ids);
+  }
+
+  friend bool operator==(const SimulationZoneID &a, const SimulationZoneID &b)
+  {
+    return a.node_ids == b.node_ids;
+  }
+};
+
+class ModifierSimulationState {
+ private:
+  mutable bool bake_loaded_;
+
+ public:
+  ModifierSimulationCache *owner_;
+  mutable std::mutex mutex_;
+  Map<SimulationZoneID, std::unique_ptr<SimulationZoneState>> zone_states_;
+  std::optional<std::string> meta_path_;
+  std::optional<std::string> bdata_dir_;
+
+  const SimulationZoneState *get_zone_state(const SimulationZoneID &zone_id) const;
+  SimulationZoneState &get_zone_state_for_write(const SimulationZoneID &zone_id);
+  void ensure_bake_loaded() const;
+};
+
+struct ModifierSimulationStateAtFrame {
+  SubFrame frame;
+  ModifierSimulationState state;
+};
+
+enum class CacheState {
+  Valid,
+  Invalid,
+  Baked,
+};
+
+struct StatesAroundFrame {
+  const ModifierSimulationStateAtFrame *prev = nullptr;
+  const ModifierSimulationStateAtFrame *current = nullptr;
+  const ModifierSimulationStateAtFrame *next = nullptr;
+};
+
+class ModifierSimulationCache {
+ private:
+  Vector<std::unique_ptr<ModifierSimulationStateAtFrame>> states_at_frames_;
+  std::unique_ptr<BDataSharing> bdata_sharing_;
+
+  friend ModifierSimulationState;
+
+ public:
+  CacheState cache_state_ = CacheState::Valid;
+  bool failed_finding_bake_ = false;
+
+  void try_discover_bake(StringRefNull meta_dir, StringRefNull bdata_dir);
+
+  bool has_state_at_frame(const SubFrame &frame) const;
+  bool has_states() const;
+  const ModifierSimulationState *get_state_at_exact_frame(const SubFrame &frame) const;
+  ModifierSimulationState &get_state_at_frame_for_write(const SubFrame &frame);
+  StatesAroundFrame get_states_around_frame(const SubFrame &frame) const;
+
+  void invalidate()
+  {
+    cache_state_ = CacheState::Invalid;
+  }
+
+  CacheState cache_state() const
+  {
+    return cache_state_;
+  }
+
+  void reset();
+};
+
+}  // namespace blender::bke::sim

source/blender/blenkernel/BKE_simulation_state_serialize.hh

@@ -0,0 +1,170 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "BKE_simulation_state.hh"
+
+#include "BLI_serialize.hh"
+
+struct Main;
+struct ModifierData;
+
+namespace blender {
+class fstream;
+}
+
+namespace blender::bke::sim {
+
+using DictionaryValue = io::serialize::DictionaryValue;
+using DictionaryValuePtr = std::shared_ptr<DictionaryValue>;
+
+/**
+ * Reference to a slice of memory typically stored on disk.
+ */
+struct BDataSlice {
+  std::string name;
+  IndexRange range;
+
+  DictionaryValuePtr serialize() const;
+  static std::optional<BDataSlice> deserialize(const io::serialize::DictionaryValue &io_slice);
+};
+
+/**
+ * Abstract base class for loading binary data.
+ */
+class BDataReader {
+ public:
+  /**
+   * Read the data from the given slice into the provided memory buffer.
+   * \return True on success, otherwise false.
+   */
+  [[nodiscard]] virtual bool read(const BDataSlice &slice, void *r_data) const = 0;
+};
+
+/**
+ * Abstract base class for writing binary data.
+ */
+class BDataWriter {
+ public:
+  /**
+   * Write the provided binary data.
+   * \return Slice where the data has been written to.
+   */
+  virtual BDataSlice write(const void *data, int64_t size) = 0;
+};
+
+/**
+ * Allows for simple data deduplication when writing or reading data by making use of implicit
+ * sharing.
+ */
+class BDataSharing {
+ private:
+  struct StoredByRuntimeValue {
+    /**
+     * Version of the shared data that was written before. This is needed because the data might
+     * be changed later without changing the #ImplicitSharingInfo pointer.
+     */
+    int64_t sharing_info_version;
+    /**
+     * Identifier of the stored data. This includes information for where the data is stored (a
+     * #BDataSlice) and optionally information for how it is loaded (e.g. endian information).
+     */
+    DictionaryValuePtr io_data;
+  };
+
+  /**
+   * Map used to detect when some data has already been written. It keeps a weak reference to
+   * #ImplicitSharingInfo, allowing it to check for equality of two arrays just by comparing the
+   * sharing info's pointer and version.
+   */
+  Map<const ImplicitSharingInfo *, StoredByRuntimeValue> stored_by_runtime_;
+
+  /**
+   * Use a mutex so that #read_shared can be implemented in a thread-safe way.
+   */
+  mutable std::mutex mutex_;
+  /**
+   * Map used to detect when some data has been previously loaded. This keeps strong
+   * references to #ImplicitSharingInfo.
+   */
+  mutable Map<std::string, ImplicitSharingInfoAndData> runtime_by_stored_;
+
+ public:
+  ~BDataSharing();
+
+  /**
+   * Check if the data referenced by `sharing_info` has been written before. If yes, return the
+   * identifier for the previously written data. Otherwise, write the data now and store the
+   * identifier for later use.
+   * \return Identifier that indicates from where the data has been written.
+   */
+  [[nodiscard]] DictionaryValuePtr write_shared(const ImplicitSharingInfo *sharing_info,
+                                                FunctionRef<DictionaryValuePtr()> write_fn);
+
+  /**
+   * Check if the data identified by `io_data` has been read before or load it now.
+   * \return Shared ownership to the read data, or none if there was an error.
+   */
+  [[nodiscard]] std::optional<ImplicitSharingInfoAndData> read_shared(
+      const DictionaryValue &io_data,
+      FunctionRef<std::optional<ImplicitSharingInfoAndData>()> read_fn) const;
+};
+
+/**
+ * A specific #BDataReader that reads from disk.
+ */
+class DiskBDataReader : public BDataReader {
+ private:
+  const std::string bdata_dir_;
+  mutable std::mutex mutex_;
+  mutable Map<std::string, std::unique_ptr<fstream>> open_input_streams_;
+
+ public:
+  DiskBDataReader(std::string bdata_dir);
+  [[nodiscard]] bool read(const BDataSlice &slice, void *r_data) const override;
+};
+
+/**
+ * A specific #BDataWriter that writes to a file on disk.
+ */
+class DiskBDataWriter : public BDataWriter {
+ private:
+  /** Name of the file that data is written to. */
+  std::string bdata_name_;
+  /** File handle. */
+  std::ostream &bdata_file_;
+  /** Current position in the file. */
+  int64_t current_offset_;
+
+ public:
+  DiskBDataWriter(std::string bdata_name, std::ostream &bdata_file, int64_t current_offset);
+
+  BDataSlice write(const void *data, int64_t size) override;
+};
+
+/**
+ * Get the directory that contains all baked simulation data for the given modifier. This is a
+ * parent directory of the two directories below.
+ */
+std::string get_bake_directory(const Main &bmain, const Object &object, const ModifierData &md);
+std::string get_bdata_directory(const Main &bmain, const Object &object, const ModifierData &md);
+std::string get_meta_directory(const Main &bmain, const Object &object, const ModifierData &md);
+
+/**
+ * Encode the simulation state in a #DictionaryValue which also contains references to external
+ * binary data that has been written using #bdata_writer.
+ */
+void serialize_modifier_simulation_state(const ModifierSimulationState &state,
+                                         BDataWriter &bdata_writer,
+                                         BDataSharing &bdata_sharing,
+                                         DictionaryValue &r_io_root);
+/**
+ * Fill the simulation state by parsing the provided #DictionaryValue which also contains
+ * references to external binary data that is read using #bdata_reader.
+ */
+void deserialize_modifier_simulation_state(const DictionaryValue &io_root,
+                                           const BDataReader &bdata_reader,
+                                           const BDataSharing &bdata_sharing,
+                                           ModifierSimulationState &r_state);
+
+}  // namespace blender::bke::sim

source/blender/blenkernel/CMakeLists.txt

@@ -232,6 +232,7 @@ set(SRC
   intern/node_tree_anonymous_attributes.cc
   intern/node_tree_field_inferencing.cc
   intern/node_tree_update.cc
+  intern/node_tree_zones.cc
   intern/object.cc
   intern/object_deform.c
   intern/object_dupli.cc
@@ -265,6 +266,8 @@ set(SRC
   intern/shader_fx.c
   intern/shrinkwrap.cc
   intern/simulation.cc
+  intern/simulation_state.cc
+  intern/simulation_state_serialize.cc
   intern/softbody.c
   intern/sound.c
   intern/speaker.c
@@ -442,6 +445,7 @@ set(SRC
   BKE_node.h
   BKE_node_runtime.hh
   BKE_node_tree_update.h
+  BKE_node_tree_zones.hh
   BKE_object.h
   BKE_object_deform.h
   BKE_object_facemap.h
@@ -464,6 +468,8 @@ set(SRC
   BKE_shader_fx.h
   BKE_shrinkwrap.h
   BKE_simulation.h
+  BKE_simulation_state.hh
+  BKE_simulation_state_serialize.hh
   BKE_softbody.h
   BKE_sound.h
   BKE_speaker.h

source/blender/blenkernel/intern/geometry_set.cc

@@ -281,6 +281,14 @@ void GeometrySet::ensure_owns_direct_data()
   }
 }
 
+void GeometrySet::ensure_owns_all_data()
+{
+  if (Instances *instances = this->get_instances_for_write()) {
+    instances->ensure_geometry_instances();
+  }
+  this->ensure_owns_direct_data();
+}
+
 bool GeometrySet::owns_direct_data() const
 {
   for (const GeometryComponentPtr &component_ptr : components_) {

source/blender/blenkernel/intern/node.cc

@@ -65,6 +65,7 @@
 #include "BKE_node.h"
 #include "BKE_node_runtime.hh"
 #include "BKE_node_tree_update.h"
+#include "BKE_node_tree_zones.hh"
 #include "BKE_type_conversions.hh"
 
 #include "RNA_access.h"
@@ -175,6 +176,10 @@ static void ntree_copy_data(Main * /*bmain*/, ID *id_dst, const ID *id_src, cons
     }
   }
 
+  for (bNode *node : ntree_dst->all_nodes()) {
+    nodeDeclarationEnsure(ntree_dst, node);
+  }
+
   /* copy interface sockets */
   BLI_listbase_clear(&ntree_dst->inputs);
   LISTBASE_FOREACH (const bNodeSocket *, src_socket, &ntree_src->inputs) {
@@ -609,6 +614,14 @@ void ntreeBlendWrite(BlendWriter *writer, bNodeTree *ntree)
         BLO_write_struct(writer, NodeImageLayer, sock->storage);
       }
     }
+    if (node->type == GEO_NODE_SIMULATION_OUTPUT) {
+      const NodeGeometrySimulationOutput &storage =
+          *static_cast<const NodeGeometrySimulationOutput *>(node->storage);
+      BLO_write_struct_array(writer, NodeSimulationItem, storage.items_num, storage.items);
+      for (const NodeSimulationItem &item : Span(storage.items, storage.items_num)) {
+        BLO_write_string(writer, item.name);
+      }
+    }
   }
 
   LISTBASE_FOREACH (bNodeLink *, link, &ntree->links) {
@@ -784,6 +797,16 @@ void ntreeBlendReadData(BlendDataReader *reader, ID *owner_id, bNodeTree *ntree)
           BLO_read_data_address(reader, &storage->string);
           break;
         }
+        case GEO_NODE_SIMULATION_OUTPUT: {
+          NodeGeometrySimulationOutput &storage = *static_cast<NodeGeometrySimulationOutput *>(
+              node->storage);
+          BLO_read_data_address(reader, &storage.items);
+          for (const NodeSimulationItem &item : Span(storage.items, storage.items_num)) {
+            BLO_read_data_address(reader, &item.name);
+          }
+          break;
+        }
+
         default:
           break;
       }
@@ -1139,7 +1162,12 @@ static void node_init(const bContext *C, bNodeTree *ntree, bNode *node)
   BLI_strncpy(node->name, DATA_(ntype->ui_name), NODE_MAXSTR);
   nodeUniqueName(ntree, node);
 
-  node_add_sockets_from_type(ntree, node, ntype);
+  /* Generally sockets should be added after the initialization, because the set of sockets might
+   * depend on node properties. */
+  const bool add_sockets_before_init = node->type == CMP_NODE_R_LAYERS;
+  if (add_sockets_before_init) {
+    node_add_sockets_from_type(ntree, node, ntype);
+  }
 
   if (ntype->initfunc != nullptr) {
     ntype->initfunc(ntree, node);
@@ -1149,6 +1177,10 @@ static void node_init(const bContext *C, bNodeTree *ntree, bNode *node)
     ntree->typeinfo->node_add_init(ntree, node);
   }
 
+  if (!add_sockets_before_init) {
+    node_add_sockets_from_type(ntree, node, ntype);
+  }
+
   if (node->id) {
     id_us_plus(node->id);
   }
@@ -2300,7 +2332,8 @@ bNode *nodeAddNode(const bContext *C, bNodeTree *ntree, const char *idname)
 
   BKE_ntree_update_tag_node_new(ntree, node);
 
-  if (ELEM(node->type, GEO_NODE_INPUT_SCENE_TIME, GEO_NODE_SELF_OBJECT)) {
+  if (ELEM(node->type, GEO_NODE_INPUT_SCENE_TIME, GEO_NODE_SELF_OBJECT, GEO_NODE_SIMULATION_INPUT))
+  {
     DEG_relations_tag_update(CTX_data_main(C));
   }
 
@@ -2424,11 +2457,6 @@ bNode *node_copy_with_mapping(bNodeTree *dst_tree,
     node_dst->typeinfo->copyfunc_api(&ptr, &node_src);
   }
 
-  /* Reset the declaration of the new node in real tree. */
-  if (dst_tree != nullptr) {
-    nodeDeclarationEnsure(dst_tree, node_dst);
-  }
-
   return node_dst;
 }
 
@@ -3265,7 +3293,9 @@ void nodeRemoveNode(Main *bmain, bNodeTree *ntree, bNode *node, const bool do_id
 
   /* Also update relations for the scene time node, which causes a dependency
    * on time that users expect to be removed when the node is removed. */
-  if (node_has_id || ELEM(node->type, GEO_NODE_INPUT_SCENE_TIME, GEO_NODE_SELF_OBJECT)) {
+  if (node_has_id ||
+      ELEM(node->type, GEO_NODE_INPUT_SCENE_TIME, GEO_NODE_SELF_OBJECT, GEO_NODE_SIMULATION_INPUT))
+  {
     if (bmain != nullptr) {
       DEG_relations_tag_update(bmain);
     }
@@ -3826,8 +3856,7 @@ bool nodeDeclarationEnsureOnOutdatedNode(bNodeTree *ntree, bNode *node)
   if (node->typeinfo->declare_dynamic) {
     BLI_assert(ntree != nullptr);
     BLI_assert(node != nullptr);
-    node->runtime->declaration = new blender::nodes::NodeDeclaration();
-    blender::nodes::build_node_declaration_dynamic(*ntree, *node, *node->runtime->declaration);
+    blender::nodes::update_node_declaration_and_sockets(*ntree, *node);
     return true;
   }
   if (node->typeinfo->declare) {

source/blender/blenkernel/intern/node_runtime.cc

@@ -255,7 +255,37 @@ struct ToposortNodeState {
   bool is_in_stack = false;
 };
 
-static void toposort_from_start_node(const ToposortDirection direction,
+static Vector<const bNode *> get_implicit_origin_nodes(const bNodeTree &ntree, bNode &node)
+{
+  Vector<const bNode *> origin_nodes;
+  if (node.type == GEO_NODE_SIMULATION_OUTPUT) {
+    for (const bNode *sim_input_node :
+         ntree.runtime->nodes_by_type.lookup(nodeTypeFind("GeometryNodeSimulationInput")))
+    {
+      const auto &storage = *static_cast<const NodeGeometrySimulationInput *>(
+          sim_input_node->storage);
+      if (storage.output_node_id == node.identifier) {
+        origin_nodes.append(sim_input_node);
+      }
+    }
+  }
+  return origin_nodes;
+}
+
+static Vector<const bNode *> get_implicit_target_nodes(const bNodeTree &ntree, bNode &node)
+{
+  Vector<const bNode *> target_nodes;
+  if (node.type == GEO_NODE_SIMULATION_INPUT) {
+    const auto &storage = *static_cast<const NodeGeometrySimulationInput *>(node.storage);
+    if (const bNode *sim_output_node = ntree.node_by_id(storage.output_node_id)) {
+      target_nodes.append(sim_output_node);
+    }
+  }
+  return target_nodes;
+}
+
+static void toposort_from_start_node(const bNodeTree &ntree,
+                                     const ToposortDirection direction,
                                      bNode &start_node,
                                      MutableSpan<ToposortNodeState> node_states,
                                      Vector<bNode *> &r_sorted_nodes,
@@ -265,6 +295,7 @@ static void toposort_from_start_node(const ToposortDirection direction,
     bNode *node;
     int socket_index = 0;
     int link_index = 0;
+    int implicit_link_index = 0;
   };
 
   Stack<Item, 64> nodes_to_check;
@@ -273,6 +304,25 @@ static void toposort_from_start_node(const ToposortDirection direction,
   while (!nodes_to_check.is_empty()) {
     Item &item = nodes_to_check.peek();
     bNode &node = *item.node;
+    bool pushed_node = false;
+
+    auto handle_linked_node = [&](bNode &linked_node) {
+      ToposortNodeState &linked_node_state = node_states[linked_node.index()];
+      if (linked_node_state.is_done) {
+        /* The linked node has already been visited. */
+        return true;
+      }
+      if (linked_node_state.is_in_stack) {
+        r_cycle_detected = true;
+      }
+      else {
+        nodes_to_check.push({&linked_node});
+        linked_node_state.is_in_stack = true;
+        pushed_node = true;
+      }
+      return false;
+    };
+
     const Span<bNodeSocket *> sockets = (direction == ToposortDirection::LeftToRight) ?
                                             node.runtime->inputs :
                                             node.runtime->outputs;
@@ -297,24 +347,38 @@ static void toposort_from_start_node(const ToposortDirection direction,
       }
       bNodeSocket &linked_socket = *socket.runtime->directly_linked_sockets[item.link_index];
       bNode &linked_node = *linked_socket.runtime->owner_node;
-      ToposortNodeState &linked_node_state = node_states[linked_node.index()];
-      if (linked_node_state.is_done) {
+      if (handle_linked_node(linked_node)) {
         /* The linked node has already been visited. */
         item.link_index++;
         continue;
       }
-      if (linked_node_state.is_in_stack) {
-        r_cycle_detected = true;
-      }
-      else {
-        nodes_to_check.push({&linked_node});
-        linked_node_state.is_in_stack = true;
-      }
       break;
     }
 
-    /* If no other element has been pushed, the current node can be pushed to the sorted list. */
-    if (&item == &nodes_to_check.peek()) {
+    if (!pushed_node) {
+      /* Some nodes are internally linked without an explicit `bNodeLink`. The toposort should
+       * still order them correctly and find cycles. */
+      const Vector<const bNode *> implicitly_linked_nodes =
+          (direction == ToposortDirection::LeftToRight) ? get_implicit_origin_nodes(ntree, node) :
+                                                          get_implicit_target_nodes(ntree, node);
+      while (true) {
+        if (item.implicit_link_index == implicitly_linked_nodes.size()) {
+          /* All implicitly linked nodes have already been visited. */
+          break;
+        }
+        const bNode &linked_node = *implicitly_linked_nodes[item.implicit_link_index];
+        if (handle_linked_node(const_cast<bNode &>(linked_node))) {
+          /* The implicitly linked node has already been visited. */
+          item.implicit_link_index++;
+          continue;
+        }
+        break;
+      }
+    }
+
+    /* If no other element has been pushed, the current node can be pushed to the sorted list.
+     */
+    if (!pushed_node) {
       ToposortNodeState &node_state = node_states[node.index()];
       node_state.is_done = true;
       node_state.is_in_stack = false;
@@ -347,7 +411,8 @@ static void update_toposort(const bNodeTree &ntree,
       /* Ignore non-start nodes. */
       continue;
     }
-    toposort_from_start_node(direction, *node, node_states, r_sorted_nodes, r_cycle_detected);
+    toposort_from_start_node(
+        ntree, direction, *node, node_states, r_sorted_nodes, r_cycle_detected);
   }
 
   if (r_sorted_nodes.size() < tree_runtime.nodes_by_id.size()) {
@@ -358,7 +423,8 @@ static void update_toposort(const bNodeTree &ntree,
         continue;
       }
       /* Start toposort at this node which is somewhere in the middle of a loop. */
-      toposort_from_start_node(direction, *node, node_states, r_sorted_nodes, r_cycle_detected);
+      toposort_from_start_node(
+          ntree, direction, *node, node_states, r_sorted_nodes, r_cycle_detected);
     }
   }
 
@@ -426,10 +492,10 @@ static void ensure_topology_cache(const bNodeTree &ntree)
     update_socket_vectors_and_owner_node(ntree);
     update_internal_link_inputs(ntree);
     update_directly_linked_links_and_sockets(ntree);
+    update_nodes_by_type(ntree);
     threading::parallel_invoke(
         tree_runtime.nodes_by_id.size() > 32,
         [&]() { update_logical_origins(ntree); },
-        [&]() { update_nodes_by_type(ntree); },
         [&]() { update_sockets_by_identifier(ntree); },
         [&]() {
           update_toposort(ntree,

source/blender/blenkernel/intern/node_tree_anonymous_attributes.cc

@@ -51,6 +51,44 @@ static const aal::RelationsInNode &get_relations_in_node(const bNode &node, Reso
       return geometry_relations;
     }
   }
+  if (ELEM(node.type, GEO_NODE_SIMULATION_INPUT, GEO_NODE_SIMULATION_OUTPUT)) {
+    aal::RelationsInNode &relations = scope.construct<aal::RelationsInNode>();
+    {
+      /* Add eval relations. */
+      int last_geometry_index = -1;
+      for (const int i : node.input_sockets().index_range()) {
+        const bNodeSocket &socket = node.input_socket(i);
+        if (socket.type == SOCK_GEOMETRY) {
+          last_geometry_index = i;
+        }
+        else if (socket_is_field(socket)) {
+          if (last_geometry_index != -1) {
+            relations.eval_relations.append({i, last_geometry_index});
+          }
+        }
+      }
+    }
+
+    {
+      /* Add available relations. */
+      int last_geometry_index = -1;
+      for (const int i : node.output_sockets().index_range()) {
+        const bNodeSocket &socket = node.output_socket(i);
+        if (socket.type == SOCK_GEOMETRY) {
+          last_geometry_index = i;
+        }
+        else if (socket_is_field(socket)) {
+          if (last_geometry_index == -1) {
+            relations.available_on_none.append(i);
+          }
+          else {
+            relations.available_relations.append({i, last_geometry_index});
+          }
+        }
+      }
+    }
+    return relations;
+  }
   if (const NodeDeclaration *node_decl = node.declaration()) {
     if (const aal::RelationsInNode *relations = node_decl->anonymous_attribute_relations()) {
       return *relations;

source/blender/blenkernel/intern/node_tree_field_inferencing.cc

@@ -258,6 +258,103 @@ static OutputFieldDependency find_group_output_dependencies(
   return OutputFieldDependency::ForPartiallyDependentField(std::move(linked_input_indices));
 }
 
+/** Result of syncing two field states. */
+enum class eFieldStateSyncResult : char {
+  /* Nothing changed. */
+  NONE = 0,
+  /* State A has been modified. */
+  CHANGED_A = (1 << 0),
+  /* State B has been modified. */
+  CHANGED_B = (1 << 1),
+};
+ENUM_OPERATORS(eFieldStateSyncResult, eFieldStateSyncResult::CHANGED_B)
+
+/**
+ * Compare both field states and select the most compatible.
+ * Afterwards both field states will be the same.
+ * \return eFieldStateSyncResult flags indicating which field states have changed.
+ */
+static eFieldStateSyncResult sync_field_states(SocketFieldState &a, SocketFieldState &b)
+{
+  const bool requires_single = a.requires_single || b.requires_single;
+  const bool is_single = a.is_single && b.is_single;
+
+  eFieldStateSyncResult res = eFieldStateSyncResult::NONE;
+  if (a.requires_single != requires_single || a.is_single != is_single) {
+    res |= eFieldStateSyncResult::CHANGED_A;
+  }
+  if (b.requires_single != requires_single || b.is_single != is_single) {
+    res |= eFieldStateSyncResult::CHANGED_B;
+  }
+
+  a.requires_single = requires_single;
+  b.requires_single = requires_single;
+  a.is_single = is_single;
+  b.is_single = is_single;
+
+  return res;
+}
+
+/**
+ * Compare field states of simulation nodes sockets and select the most compatible.
+ * Afterwards all field states will be the same.
+ * \return eFieldStateSyncResult flags indicating which field states have changed.
+ */
+static eFieldStateSyncResult simulation_nodes_field_state_sync(
+    const bNode &input_node,
+    const bNode &output_node,
+    const MutableSpan<SocketFieldState> field_state_by_socket_id)
+{
+  eFieldStateSyncResult res = eFieldStateSyncResult::NONE;
+  for (const int i : output_node.output_sockets().index_range()) {
+    /* First input node output is Delta Time which does not appear in the output node outputs. */
+    const bNodeSocket &input_socket = input_node.output_socket(i + 1);
+    const bNodeSocket &output_socket = output_node.output_socket(i);
+    SocketFieldState &input_state = field_state_by_socket_id[input_socket.index_in_tree()];
+    SocketFieldState &output_state = field_state_by_socket_id[output_socket.index_in_tree()];
+    res |= sync_field_states(input_state, output_state);
+  }
+  return res;
+}
+
+static bool propagate_special_data_requirements(
+    const bNodeTree &tree,
+    const bNode &node,
+    const MutableSpan<SocketFieldState> field_state_by_socket_id)
+{
+  tree.ensure_topology_cache();
+
+  bool need_update = false;
+
+  /* Sync field state between simulation nodes and schedule another pass if necessary. */
+  if (node.type == GEO_NODE_SIMULATION_INPUT) {
+    const NodeGeometrySimulationInput &data = *static_cast<const NodeGeometrySimulationInput *>(
+        node.storage);
+    if (const bNode *output_node = tree.node_by_id(data.output_node_id)) {
+      const eFieldStateSyncResult sync_result = simulation_nodes_field_state_sync(
+          node, *output_node, field_state_by_socket_id);
+      if (bool(sync_result & eFieldStateSyncResult::CHANGED_B)) {
+        need_update = true;
+      }
+    }
+  }
+  else if (node.type == GEO_NODE_SIMULATION_OUTPUT) {
+    for (const bNode *input_node : tree.nodes_by_type("GeometryNodeSimulationInput")) {
+      const NodeGeometrySimulationInput &data = *static_cast<const NodeGeometrySimulationInput *>(
+          input_node->storage);
+      if (node.identifier == data.output_node_id) {
+        const eFieldStateSyncResult sync_result = simulation_nodes_field_state_sync(
+            *input_node, node, field_state_by_socket_id);
+        if (bool(sync_result & eFieldStateSyncResult::CHANGED_A)) {
+          need_update = true;
+        }
+      }
+    }
+  }
+
+  return need_update;
+}
+
 static void propagate_data_requirements_from_right_to_left(
     const bNodeTree &tree,
     const Span<const FieldInferencingInterface *> interface_by_node,
@@ -265,70 +362,85 @@ static void propagate_data_requirements_from_right_to_left(
 {
   const Span<const bNode *> toposort_result = tree.toposort_right_to_left();
 
-  for (const bNode *node : toposort_result) {
-    const FieldInferencingInterface &inferencing_interface = *interface_by_node[node->index()];
+  while (true) {
+    /* Node updates may require sevaral passes due to cyclic dependencies caused by simulation
+     * input/output nodes. */
+    bool need_update = false;
 
-    for (const bNodeSocket *output_socket : node->output_sockets()) {
-      SocketFieldState &state = field_state_by_socket_id[output_socket->index_in_tree()];
+    for (const bNode *node : toposort_result) {
+      const FieldInferencingInterface &inferencing_interface = *interface_by_node[node->index()];
 
-      const OutputFieldDependency &field_dependency =
-          inferencing_interface.outputs[output_socket->index()];
+      for (const bNodeSocket *output_socket : node->output_sockets()) {
+        SocketFieldState &state = field_state_by_socket_id[output_socket->index_in_tree()];
 
-      if (field_dependency.field_type() == OutputSocketFieldType::FieldSource) {
-        continue;
-      }
-      if (field_dependency.field_type() == OutputSocketFieldType::None) {
-        state.requires_single = true;
-        state.is_always_single = true;
-        continue;
-      }
+        const OutputFieldDependency &field_dependency =
+            inferencing_interface.outputs[output_socket->index()];
 
-      /* The output is required to be a single value when it is connected to any input that does
-       * not support fields. */
-      for (const bNodeSocket *target_socket : output_socket->directly_linked_sockets()) {
-        if (target_socket->is_available()) {
-          state.requires_single |=
-              field_state_by_socket_id[target_socket->index_in_tree()].requires_single;
+        if (field_dependency.field_type() == OutputSocketFieldType::FieldSource) {
+          continue;
+        }
+        if (field_dependency.field_type() == OutputSocketFieldType::None) {
+          state.requires_single = true;
+          state.is_always_single = true;
+          continue;
         }
-      }
 
-      if (state.requires_single) {
-        bool any_input_is_field_implicitly = false;
-        const Vector<const bNodeSocket *> connected_inputs = gather_input_socket_dependencies(
-            field_dependency, *node);
-        for (const bNodeSocket *input_socket : connected_inputs) {
-          if (!input_socket->is_available()) {
-            continue;
+        /* The output is required to be a single value when it is connected to any input that does
+         * not support fields. */
+        for (const bNodeSocket *target_socket : output_socket->directly_linked_sockets()) {
+          if (target_socket->is_available()) {
+            state.requires_single |=
+                field_state_by_socket_id[target_socket->index_in_tree()].requires_single;
           }
-          if (inferencing_interface.inputs[input_socket->index()] ==
-              InputSocketFieldType::Implicit) {
-            if (!input_socket->is_logically_linked()) {
-              any_input_is_field_implicitly = true;
-              break;
+        }
+
+        if (state.requires_single) {
+          bool any_input_is_field_implicitly = false;
+          const Vector<const bNodeSocket *> connected_inputs = gather_input_socket_dependencies(
+              field_dependency, *node);
+          for (const bNodeSocket *input_socket : connected_inputs) {
+            if (!input_socket->is_available()) {
+              continue;
+            }
+            if (inferencing_interface.inputs[input_socket->index()] ==
+                InputSocketFieldType::Implicit) {
+              if (!input_socket->is_logically_linked()) {
+                any_input_is_field_implicitly = true;
+                break;
+              }
+            }
+          }
+          if (any_input_is_field_implicitly) {
+            /* This output isn't a single value actually. */
+            state.requires_single = false;
+          }
+          else {
+            /* If the output is required to be a single value, the connected inputs in the same
+             * node must not be fields as well. */
+            for (const bNodeSocket *input_socket : connected_inputs) {
+              field_state_by_socket_id[input_socket->index_in_tree()].requires_single = true;
             }
           }
         }
-        if (any_input_is_field_implicitly) {
-          /* This output isn't a single value actually. */
-          state.requires_single = false;
-        }
-        else {
-          /* If the output is required to be a single value, the connected inputs in the same node
-           * must not be fields as well. */
-          for (const bNodeSocket *input_socket : connected_inputs) {
-            field_state_by_socket_id[input_socket->index_in_tree()].requires_single = true;
-          }
+      }
+
+      /* Some inputs do not require fields independent of what the outputs are connected to. */
+      for (const bNodeSocket *input_socket : node->input_sockets()) {
+        SocketFieldState &state = field_state_by_socket_id[input_socket->index_in_tree()];
+        if (inferencing_interface.inputs[input_socket->index()] == InputSocketFieldType::None) {
+          state.requires_single = true;
+          state.is_always_single = true;
         }
       }
+
+      /* Find reverse dependencies and resolve conflicts, which may require another pass. */
+      if (propagate_special_data_requirements(tree, *node, field_state_by_socket_id)) {
+        need_update = true;
+      }
     }
 
-    /* Some inputs do not require fields independent of what the outputs are connected to. */
-    for (const bNodeSocket *input_socket : node->input_sockets()) {
-      SocketFieldState &state = field_state_by_socket_id[input_socket->index_in_tree()];
-      if (inferencing_interface.inputs[input_socket->index()] == InputSocketFieldType::None) {
-        state.requires_single = true;
-        state.is_always_single = true;
-      }
+    if (!need_update) {
+      break;
     }
   }
 }
@@ -384,68 +496,82 @@ static void propagate_field_status_from_left_to_right(
 {
   const Span<const bNode *> toposort_result = tree.toposort_left_to_right();
 
-  for (const bNode *node : toposort_result) {
-    if (node->type == NODE_GROUP_INPUT) {
-      continue;
-    }
+  while (true) {
+    /* Node updates may require sevaral passes due to cyclic dependencies. */
+    bool need_update = false;
 
-    const FieldInferencingInterface &inferencing_interface = *interface_by_node[node->index()];
-
-    /* Update field state of input sockets, also taking into account linked origin sockets. */
-    for (const bNodeSocket *input_socket : node->input_sockets()) {
-      SocketFieldState &state = field_state_by_socket_id[input_socket->index_in_tree()];
-      if (state.is_always_single) {
-        state.is_single = true;
+    for (const bNode *node : toposort_result) {
+      if (node->type == NODE_GROUP_INPUT) {
         continue;
       }
-      state.is_single = true;
-      if (!input_socket->is_directly_linked()) {
-        if (inferencing_interface.inputs[input_socket->index()] == InputSocketFieldType::Implicit)
-        {
-          state.is_single = false;
+
+      const FieldInferencingInterface &inferencing_interface = *interface_by_node[node->index()];
+
+      /* Update field state of input sockets, also taking into account linked origin sockets. */
+      for (const bNodeSocket *input_socket : node->input_sockets()) {
+        SocketFieldState &state = field_state_by_socket_id[input_socket->index_in_tree()];
+        if (state.is_always_single) {
+          state.is_single = true;
+          continue;
         }
-      }
-      else {
-        for (const bNodeSocket *origin_socket : input_socket->directly_linked_sockets()) {
-          if (!field_state_by_socket_id[origin_socket->index_in_tree()].is_single) {
+        state.is_single = true;
+        if (!input_socket->is_directly_linked()) {
+          if (inferencing_interface.inputs[input_socket->index()] ==
+              InputSocketFieldType::Implicit) {
             state.is_single = false;
-            break;
           }
         }
-      }
-    }
-
-    /* Update field state of output sockets, also taking into account input sockets. */
-    for (const bNodeSocket *output_socket : node->output_sockets()) {
-      SocketFieldState &state = field_state_by_socket_id[output_socket->index_in_tree()];
-      const OutputFieldDependency &field_dependency =
-          inferencing_interface.outputs[output_socket->index()];
-
-      switch (field_dependency.field_type()) {
-        case OutputSocketFieldType::None: {
-          state.is_single = true;
-          break;
-        }
-        case OutputSocketFieldType::FieldSource: {
-          state.is_single = false;
-          state.is_field_source = true;
-          break;
-        }
-        case OutputSocketFieldType::PartiallyDependent:
-        case OutputSocketFieldType::DependentField: {
-          for (const bNodeSocket *input_socket :
-               gather_input_socket_dependencies(field_dependency, *node)) {
-            if (!input_socket->is_available()) {
-              continue;
-            }
-            if (!field_state_by_socket_id[input_socket->index_in_tree()].is_single) {
+        else {
+          for (const bNodeSocket *origin_socket : input_socket->directly_linked_sockets()) {
+            if (!field_state_by_socket_id[origin_socket->index_in_tree()].is_single) {
               state.is_single = false;
               break;
             }
           }
-          break;
         }
       }
+
+      /* Update field state of output sockets, also taking into account input sockets. */
+      for (const bNodeSocket *output_socket : node->output_sockets()) {
+        SocketFieldState &state = field_state_by_socket_id[output_socket->index_in_tree()];
+        const OutputFieldDependency &field_dependency =
+            inferencing_interface.outputs[output_socket->index()];
+
+        switch (field_dependency.field_type()) {
+          case OutputSocketFieldType::None: {
+            state.is_single = true;
+            break;
+          }
+          case OutputSocketFieldType::FieldSource: {
+            state.is_single = false;
+            state.is_field_source = true;
+            break;
+          }
+          case OutputSocketFieldType::PartiallyDependent:
+          case OutputSocketFieldType::DependentField: {
+            for (const bNodeSocket *input_socket :
+                 gather_input_socket_dependencies(field_dependency, *node)) {
+              if (!input_socket->is_available()) {
+                continue;
+              }
+              if (!field_state_by_socket_id[input_socket->index_in_tree()].is_single) {
+                state.is_single = false;
+                break;
+              }
+            }
+            break;
+          }
+        }
+      }
+
+      /* Find reverse dependencies and resolve conflicts, which may require another pass. */
+      if (propagate_special_data_requirements(tree, *node, field_state_by_socket_id)) {
+        need_update = true;
+      }
+    }
+
+    if (!need_update) {
+      break;
     }
   }
 }

source/blender/blenkernel/intern/node_tree_update.cc

@@ -53,6 +53,7 @@ static void add_tree_tag(bNodeTree *ntree, const eNodeTreeChangedFlag flag)
 {
   ntree->runtime->changed_flag |= flag;
   ntree->runtime->topology_cache_mutex.tag_dirty();
+  ntree->runtime->tree_zones_cache_mutex.tag_dirty();
 }
 
 static void add_node_tag(bNodeTree *ntree, bNode *node, const eNodeTreeChangedFlag flag)
@@ -572,6 +573,16 @@ class NodeTreeMainUpdater {
         return true;
       }
     }
+    /* Check paired simulation zone nodes. */
+    if (node.type == GEO_NODE_SIMULATION_INPUT) {
+      const NodeGeometrySimulationInput *data = static_cast<const NodeGeometrySimulationInput *>(
+          node.storage);
+      if (const bNode *output_node = ntree.node_by_id(data->output_node_id)) {
+        if (output_node->runtime->changed_flag & NTREE_CHANGED_NODE_PROPERTY) {
+          return true;
+        }
+      }
+    }
     return false;
   }
 

source/blender/blenkernel/intern/node_tree_zones.cc

@@ -0,0 +1,252 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BKE_node.h"
+#include "BKE_node_runtime.hh"
+#include "BKE_node_tree_zones.hh"
+
+#include "BLI_bit_group_vector.hh"
+#include "BLI_bit_span_ops.hh"
+#include "BLI_task.hh"
+#include "BLI_timeit.hh"
+
+namespace blender::bke::node_tree_zones {
+
+static void update_zone_depths(TreeZone &zone)
+{
+  if (zone.depth >= 0) {
+    return;
+  }
+  if (zone.parent_zone == nullptr) {
+    zone.depth = 0;
+    return;
+  }
+  update_zone_depths(*zone.parent_zone);
+  zone.depth = zone.parent_zone->depth + 1;
+}
+
+static Vector<std::unique_ptr<TreeZone>> find_zone_nodes(
+    const bNodeTree &tree, TreeZones &owner, Map<const bNode *, TreeZone *> &r_zone_by_inout_node)
+{
+  Vector<std::unique_ptr<TreeZone>> zones;
+  for (const bNode *node : tree.nodes_by_type("GeometryNodeSimulationOutput")) {
+    auto zone = std::make_unique<TreeZone>();
+    zone->owner = &owner;
+    zone->index = zones.size();
+    zone->output_node = node;
+    r_zone_by_inout_node.add(node, zone.get());
+    zones.append_and_get_index(std::move(zone));
+  }
+  for (const bNode *node : tree.nodes_by_type("GeometryNodeSimulationInput")) {
+    const auto &storage = *static_cast<NodeGeometrySimulationInput *>(node->storage);
+    if (const bNode *sim_output_node = tree.node_by_id(storage.output_node_id)) {
+      if (TreeZone *zone = r_zone_by_inout_node.lookup_default(sim_output_node, nullptr)) {
+        zone->input_node = node;
+        r_zone_by_inout_node.add(node, zone);
+      }
+    }
+  }
+  return zones;
+}
+
+struct ZoneRelation {
+  TreeZone *parent;
+  TreeZone *child;
+};
+
+static Vector<ZoneRelation> get_direct_zone_relations(
+    const Span<std::unique_ptr<TreeZone>> all_zones,
+    const BitGroupVector<> &depend_on_input_flag_array)
+{
+  Vector<ZoneRelation> zone_relations;
+
+  /* Gather all relations, even the transitive once. */
+  for (const std::unique_ptr<TreeZone> &zone : all_zones) {
+    const int zone_i = zone->index;
+    for (const bNode *node : {zone->output_node}) {
+      if (node == nullptr) {
+        continue;
+      }
+      const BoundedBitSpan depend_on_input_flags = depend_on_input_flag_array[node->index()];
+      bits::foreach_1_index(depend_on_input_flags, [&](const int parent_zone_i) {
+        if (parent_zone_i != zone_i) {
+          zone_relations.append({all_zones[parent_zone_i].get(), zone.get()});
+        }
+      });
+    }
+  }
+
+  /* Remove transitive relations. This is a brute force algorithm currently. */
+  Vector<int> transitive_relations;
+  for (const int a : zone_relations.index_range()) {
+    const ZoneRelation &relation_a = zone_relations[a];
+    for (const int b : zone_relations.index_range()) {
+      if (a == b) {
+        continue;
+      }
+      const ZoneRelation &relation_b = zone_relations[b];
+      for (const int c : zone_relations.index_range()) {
+        if (a == c || b == c) {
+          continue;
+        }
+        const ZoneRelation &relation_c = zone_relations[c];
+        if (relation_a.child == relation_b.parent && relation_a.parent == relation_c.parent &&
+            relation_b.child == relation_c.child)
+        {
+          transitive_relations.append_non_duplicates(c);
+        }
+      }
+    }
+  }
+  std::sort(transitive_relations.begin(), transitive_relations.end(), std::greater<>());
+  for (const int i : transitive_relations) {
+    zone_relations.remove_and_reorder(i);
+  }
+
+  return zone_relations;
+}
+
+static void update_parent_zone_per_node(const Span<const bNode *> all_nodes,
+                                        const Span<std::unique_ptr<TreeZone>> all_zones,
+                                        const BitGroupVector<> &depend_on_input_flag_array,
+                                        Map<int, int> &r_parent_zone_by_node_id)
+{
+  for (const int node_i : all_nodes.index_range()) {
+    const bNode &node = *all_nodes[node_i];
+    const BoundedBitSpan depend_on_input_flags = depend_on_input_flag_array[node_i];
+    TreeZone *parent_zone = nullptr;
+    bits::foreach_1_index(depend_on_input_flags, [&](const int parent_zone_i) {
+      TreeZone *zone = all_zones[parent_zone_i].get();
+      if (ELEM(&node, zone->input_node, zone->output_node)) {
+        return;
+      }
+      if (parent_zone == nullptr || zone->depth > parent_zone->depth) {
+        parent_zone = zone;
+      }
+    });
+    if (parent_zone != nullptr) {
+      r_parent_zone_by_node_id.add(node.identifier, parent_zone->index);
+    }
+  }
+}
+
+static std::unique_ptr<TreeZones> discover_tree_zones(const bNodeTree &tree)
+{
+  if (tree.has_available_link_cycle()) {
+    return {};
+  }
+
+  std::unique_ptr<TreeZones> tree_zones = std::make_unique<TreeZones>();
+
+  const Span<const bNode *> all_nodes = tree.all_nodes();
+  Map<const bNode *, TreeZone *> zone_by_inout_node;
+  tree_zones->zones = find_zone_nodes(tree, *tree_zones, zone_by_inout_node);
+
+  const int zones_num = tree_zones->zones.size();
+  const int nodes_num = all_nodes.size();
+  /* A bit for every node-zone-combination. The bit is set when the node is in the zone. */
+  BitGroupVector<> depend_on_input_flag_array(nodes_num, zones_num, false);
+  /* The bit is set when the node depends on the output of the zone. */
+  BitGroupVector<> depend_on_output_flag_array(nodes_num, zones_num, false);
+
+  const Span<const bNode *> sorted_nodes = tree.toposort_left_to_right();
+  for (const bNode *node : sorted_nodes) {
+    const int node_i = node->index();
+    MutableBoundedBitSpan depend_on_input_flags = depend_on_input_flag_array[node_i];
+    MutableBoundedBitSpan depend_on_output_flags = depend_on_output_flag_array[node_i];
+
+    /* Forward all bits from the nodes to the left. */
+    for (const bNodeSocket *input_socket : node->input_sockets()) {
+      if (!input_socket->is_available()) {
+        continue;
+      }
+      for (const bNodeLink *link : input_socket->directly_linked_links()) {
+        if (link->is_muted()) {
+          continue;
+        }
+        const bNode &from_node = *link->fromnode;
+        const int from_node_i = from_node.index();
+        depend_on_input_flags |= depend_on_input_flag_array[from_node_i];
+        depend_on_output_flags |= depend_on_output_flag_array[from_node_i];
+      }
+    }
+    if (node->type == GEO_NODE_SIMULATION_INPUT) {
+      if (const TreeZone *zone = zone_by_inout_node.lookup_default(node, nullptr)) {
+        /* Now entering a zone, so set the corresponding bit. */
+        depend_on_input_flags[zone->index].set();
+      }
+    }
+    else if (node->type == GEO_NODE_SIMULATION_OUTPUT) {
+      if (const TreeZone *zone = zone_by_inout_node.lookup_default(node, nullptr)) {
+        /* The output is implicitly linked to the input, so also propagate the bits from there. */
+        if (const bNode *zone_input_node = zone->input_node) {
+          const int input_node_i = zone_input_node->index();
+          depend_on_input_flags |= depend_on_input_flag_array[input_node_i];
+          depend_on_output_flags |= depend_on_output_flag_array[input_node_i];
+        }
+        /* Now exiting a zone, so change the bits accordingly. */
+        depend_on_input_flags[zone->index].reset();
+        depend_on_output_flags[zone->index].set();
+      }
+    }
+
+    if (bits::has_common_set_bits(depend_on_input_flags, depend_on_output_flags)) {
+      /* A node can not be inside and after a zone at the same time. */
+      return {};
+    }
+  }
+
+  const Vector<ZoneRelation> zone_relations = get_direct_zone_relations(
+      tree_zones->zones, depend_on_input_flag_array);
+
+  /* Set parent and child pointers in zones. */
+  for (const ZoneRelation &relation : zone_relations) {
+    relation.parent->child_zones.append(relation.child);
+    BLI_assert(relation.child->parent_zone == nullptr);
+    relation.child->parent_zone = relation.parent;
+  }
+
+  /* Update depths. */
+  for (std::unique_ptr<TreeZone> &zone : tree_zones->zones) {
+    update_zone_depths(*zone);
+  }
+
+  update_parent_zone_per_node(all_nodes,
+                              tree_zones->zones,
+                              depend_on_input_flag_array,
+                              tree_zones->parent_zone_by_node_id);
+
+  for (const int node_i : all_nodes.index_range()) {
+    const bNode *node = all_nodes[node_i];
+    const int parent_zone_i = tree_zones->parent_zone_by_node_id.lookup_default(node->identifier,
+                                                                                -1);
+    if (parent_zone_i != -1) {
+      tree_zones->zones[parent_zone_i]->child_nodes.append(node);
+    }
+  }
+
+  return tree_zones;
+}
+
+const TreeZones *get_tree_zones(const bNodeTree &tree)
+{
+  tree.runtime->tree_zones_cache_mutex.ensure(
+      [&]() { tree.runtime->tree_zones = discover_tree_zones(tree); });
+  return tree.runtime->tree_zones.get();
+}
+
+bool TreeZone::contains_node_recursively(const bNode &node) const
+{
+  const TreeZones *zones = this->owner;
+  const int parent_zone_i = zones->parent_zone_by_node_id.lookup_default(node.identifier, -1);
+  if (parent_zone_i == -1) {
+    return false;
+  }
+  for (const TreeZone *zone = zones->zones[parent_zone_i].get(); zone; zone = zone->parent_zone) {
+    if (zone == this) {
+      return true;
+    }
+  }
+  return false;
+}
+
+}  // namespace blender::bke::node_tree_zones

source/blender/blenkernel/intern/simulation_state.cc

@@ -0,0 +1,203 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BKE_curves.hh"
+#include "BKE_simulation_state.hh"
+#include "BKE_simulation_state_serialize.hh"
+
+#include "DNA_curves_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_node_types.h"
+#include "DNA_pointcloud_types.h"
+
+#include "BLI_fileops.hh"
+#include "BLI_hash_md5.h"
+#include "BLI_path_util.h"
+
+namespace blender::bke::sim {
+
+GeometrySimulationStateItem::GeometrySimulationStateItem(GeometrySet geometry)
+    : geometry_(std::move(geometry))
+{
+}
+
+PrimitiveSimulationStateItem::PrimitiveSimulationStateItem(const CPPType &type, const void *value)
+    : type_(type)
+{
+  value_ = MEM_mallocN_aligned(type.size(), type.alignment(), __func__);
+  type.copy_construct(value, value_);
+}
+
+PrimitiveSimulationStateItem::~PrimitiveSimulationStateItem()
+{
+  type_.destruct(value_);
+  MEM_freeN(value_);
+}
+
+StringSimulationStateItem::StringSimulationStateItem(std::string value) : value_(std::move(value))
+{
+}
+
+void ModifierSimulationCache::try_discover_bake(const StringRefNull meta_dir,
+                                                const StringRefNull bdata_dir)
+{
+  if (failed_finding_bake_) {
+    return;
+  }
+  if (!BLI_is_dir(meta_dir.c_str()) || !BLI_is_dir(bdata_dir.c_str())) {
+    failed_finding_bake_ = true;
+    return;
+  }
+
+  direntry *dir_entries = nullptr;
+  const int dir_entries_num = BLI_filelist_dir_contents(meta_dir.c_str(), &dir_entries);
+  BLI_SCOPED_DEFER([&]() { BLI_filelist_free(dir_entries, dir_entries_num); });
+
+  if (dir_entries_num == 0) {
+    failed_finding_bake_ = true;
+    return;
+  }
+
+  this->reset();
+
+  for (const int i : IndexRange(dir_entries_num)) {
+    const direntry &dir_entry = dir_entries[i];
+    const StringRefNull dir_entry_path = dir_entry.path;
+    if (!dir_entry_path.endswith(".json")) {
+      continue;
+    }
+    char modified_file_name[FILENAME_MAX];
+    BLI_strncpy(modified_file_name, dir_entry.relname, sizeof(modified_file_name));
+    BLI_str_replace_char(modified_file_name, '_', '.');
+
+    const SubFrame frame = std::stof(modified_file_name);
+
+    auto new_state_at_frame = std::make_unique<ModifierSimulationStateAtFrame>();
+    new_state_at_frame->frame = frame;
+    new_state_at_frame->state.bdata_dir_ = bdata_dir;
+    new_state_at_frame->state.meta_path_ = dir_entry.path;
+    new_state_at_frame->state.owner_ = this;
+    states_at_frames_.append(std::move(new_state_at_frame));
+  }
+
+  bdata_sharing_ = std::make_unique<BDataSharing>();
+
+  cache_state_ = CacheState::Baked;
+}
+
+bool ModifierSimulationCache::has_state_at_frame(const SubFrame &frame) const
+{
+  for (const auto &item : states_at_frames_) {
+    if (item->frame == frame) {
+      return true;
+    }
+  }
+  return false;
+}
+
+bool ModifierSimulationCache::has_states() const
+{
+  return !states_at_frames_.is_empty();
+}
+
+const ModifierSimulationState *ModifierSimulationCache::get_state_at_exact_frame(
+    const SubFrame &frame) const
+{
+  for (const auto &item : states_at_frames_) {
+    if (item->frame == frame) {
+      return &item->state;
+    }
+  }
+  return nullptr;
+}
+
+ModifierSimulationState &ModifierSimulationCache::get_state_at_frame_for_write(
+    const SubFrame &frame)
+{
+  for (const auto &item : states_at_frames_) {
+    if (item->frame == frame) {
+      return item->state;
+    }
+  }
+  states_at_frames_.append(std::make_unique<ModifierSimulationStateAtFrame>());
+  states_at_frames_.last()->frame = frame;
+  states_at_frames_.last()->state.owner_ = this;
+  return states_at_frames_.last()->state;
+}
+
+StatesAroundFrame ModifierSimulationCache::get_states_around_frame(const SubFrame &frame) const
+{
+  StatesAroundFrame states_around_frame;
+  for (const auto &item : states_at_frames_) {
+    if (item->frame < frame) {
+      if (states_around_frame.prev == nullptr || item->frame > states_around_frame.prev->frame) {
+        states_around_frame.prev = item.get();
+      }
+    }
+    if (item->frame == frame) {
+      if (states_around_frame.current == nullptr) {
+        states_around_frame.current = item.get();
+      }
+    }
+    if (item->frame > frame) {
+      if (states_around_frame.next == nullptr || item->frame < states_around_frame.next->frame) {
+        states_around_frame.next = item.get();
+      }
+    }
+  }
+  return states_around_frame;
+}
+
+const SimulationZoneState *ModifierSimulationState::get_zone_state(
+    const SimulationZoneID &zone_id) const
+{
+  std::lock_guard lock{mutex_};
+  if (auto *ptr = zone_states_.lookup_ptr(zone_id)) {
+    return ptr->get();
+  }
+  return nullptr;
+}
+
+SimulationZoneState &ModifierSimulationState::get_zone_state_for_write(
+    const SimulationZoneID &zone_id)
+{
+  std::lock_guard lock{mutex_};
+  return *zone_states_.lookup_or_add_cb(zone_id,
+                                        []() { return std::make_unique<SimulationZoneState>(); });
+}
+
+void ModifierSimulationState::ensure_bake_loaded() const
+{
+  std::scoped_lock lock{mutex_};
+  if (bake_loaded_) {
+    return;
+  }
+  if (!meta_path_ || !bdata_dir_) {
+    return;
+  }
+
+  const std::shared_ptr<io::serialize::Value> io_root_value = io::serialize::read_json_file(
+      *meta_path_);
+  if (!io_root_value) {
+    return;
+  }
+  const DictionaryValue *io_root = io_root_value->as_dictionary_value();
+  if (!io_root) {
+    return;
+  }
+
+  const DiskBDataReader bdata_reader{*bdata_dir_};
+  deserialize_modifier_simulation_state(*io_root,
+                                        bdata_reader,
+                                        *owner_->bdata_sharing_,
+                                        const_cast<ModifierSimulationState &>(*this));
+  bake_loaded_ = true;
+}
+
+void ModifierSimulationCache::reset()
+{
+  states_at_frames_.clear();
+  bdata_sharing_.reset();
+  cache_state_ = CacheState::Valid;
+}
+
+}  // namespace blender::bke::sim

source/blender/blenlib/BLI_implicit_sharing.hh

@@ -193,6 +193,14 @@ class ImplicitSharingMixin : public ImplicitSharingInfo {
   virtual void delete_self() = 0;
 };
 
+/**
+ * Utility that contains sharing information and the data that is shared.
+ */
+struct ImplicitSharingInfoAndData {
+  const ImplicitSharingInfo *sharing_info = nullptr;
+  const void *data = nullptr;
+};
+
 namespace implicit_sharing {
 
 namespace detail {

source/blender/blenlib/BLI_sub_frame.hh

@@ -0,0 +1,93 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "BLI_assert.h"
+#include "BLI_math_base.h"
+
+namespace blender {
+
+/**
+ * Contains an integer frame number and a subframe float in the range [0, 1).
+ */
+struct SubFrame {
+ private:
+  int frame_;
+  float subframe_;
+
+ public:
+  SubFrame(const int frame = 0, float subframe = 0.0f) : frame_(frame), subframe_(subframe)
+  {
+    BLI_assert(subframe >= 0.0f);
+    BLI_assert(subframe < 1.0f);
+  }
+
+  SubFrame(const float frame) : SubFrame(int(floorf(frame)), fractf(frame)) {}
+
+  int frame() const
+  {
+    return frame_;
+  }
+
+  float subframe() const
+  {
+    return subframe_;
+  }
+
+  explicit operator float() const
+  {
+    return float(frame_) + float(subframe_);
+  }
+
+  explicit operator double() const
+  {
+    return double(frame_) + double(subframe_);
+  }
+
+  static SubFrame min()
+  {
+    return {INT32_MIN, 0.0f};
+  }
+
+  static SubFrame max()
+  {
+    return {INT32_MAX, std::nexttowardf(1.0f, 0.0)};
+  }
+
+  friend bool operator==(const SubFrame &a, const SubFrame &b)
+  {
+    return a.frame_ == b.frame_ && a.subframe_ == b.subframe_;
+  }
+
+  friend bool operator!=(const SubFrame &a, const SubFrame &b)
+  {
+    return !(a == b);
+  }
+
+  friend bool operator<(const SubFrame &a, const SubFrame &b)
+  {
+    return a.frame_ < b.frame_ || (a.frame_ == b.frame_ && a.subframe_ < b.subframe_);
+  }
+
+  friend bool operator<=(const SubFrame &a, const SubFrame &b)
+  {
+    return a.frame_ <= b.frame_ || (a.frame_ == b.frame_ && a.subframe_ <= b.subframe_);
+  }
+
+  friend bool operator>(const SubFrame &a, const SubFrame &b)
+  {
+    return a.frame_ > b.frame_ || (a.frame_ == b.frame_ && a.subframe_ > b.subframe_);
+  }
+
+  friend bool operator>=(const SubFrame &a, const SubFrame &b)
+  {
+    return a.frame_ >= b.frame_ || (a.frame_ == b.frame_ && a.subframe_ >= b.subframe_);
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const SubFrame &a)
+  {
+    return stream << float(a);
+  }
+};
+
+}  // namespace blender

source/blender/blenlib/CMakeLists.txt

@@ -352,6 +352,7 @@ set(SRC
   BLI_string_search.h
   BLI_string_utf8.h
   BLI_string_utils.h
+  BLI_sub_frame.hh
   BLI_sys_types.h
   BLI_system.h
   BLI_task.h

source/blender/blenloader/intern/versioning_300.cc

@@ -4329,6 +4329,83 @@ void blo_do_versions_300(FileData *fd, Library * /*lib*/, Main *bmain)
 
     /* Rename Grease Pencil weight draw brush. */
     do_versions_rename_id(bmain, ID_BR, "Draw Weight", "Weight Draw");
+
+    /* Identifier generation for simulation node sockets changed.
+     * Update identifiers so links are not removed during validation.
+     * This only affects files that have been created using simulation nodes before they were first
+     * officially released. */
+    if (!DNA_struct_elem_find(fd->filesdna, "NodeSimulationItem", "int", "identifier")) {
+      static auto set_socket_identifiers =
+          [](bNode *node, const bNode *output_node, int extra_outputs) {
+            const NodeGeometrySimulationOutput *output_data =
+                static_cast<const NodeGeometrySimulationOutput *>(output_node->storage);
+            /* Includes extension socket. */
+            BLI_assert(BLI_listbase_count(&node->inputs) == output_data->items_num + 1);
+            /* Includes extension socket. */
+            BLI_assert(BLI_listbase_count(&node->outputs) ==
+                       output_data->items_num + 1 + extra_outputs);
+
+            int i;
+            LISTBASE_FOREACH_INDEX (bNodeSocket *, sock, &node->inputs, i) {
+              /* Skip extension socket. */
+              if (i >= output_data->items_num) {
+                break;
+              }
+              BLI_snprintf(sock->identifier,
+                           sizeof(sock->identifier),
+                           "Item_%d",
+                           output_data->items[i].identifier);
+            }
+            LISTBASE_FOREACH_INDEX (bNodeSocket *, sock, &node->outputs, i) {
+              const int item_i = i - extra_outputs;
+              /* Skip extra outputs. */
+              if (i < extra_outputs) {
+                continue;
+              }
+              /* Skip extension socket. */
+              if (item_i >= output_data->items_num) {
+                break;
+              }
+              BLI_snprintf(sock->identifier,
+                           sizeof(sock->identifier),
+                           "Item_%d",
+                           output_data->items[i - extra_outputs].identifier);
+            }
+          };
+
+      LISTBASE_FOREACH (bNodeTree *, ntree, &bmain->nodetrees) {
+        if (ntree->type == NTREE_GEOMETRY) {
+          /* Initialize item identifiers. */
+          LISTBASE_FOREACH (bNode *, node, &ntree->nodes) {
+            if (node->type == GEO_NODE_SIMULATION_OUTPUT) {
+              NodeGeometrySimulationOutput *data = static_cast<NodeGeometrySimulationOutput *>(
+                  node->storage);
+              for (int i = 0; i < data->items_num; ++i) {
+                data->items[i].identifier = i;
+              }
+              data->next_identifier = data->items_num;
+            }
+          }
+
+          LISTBASE_FOREACH (bNode *, node, &ntree->nodes) {
+            if (node->type == GEO_NODE_SIMULATION_INPUT) {
+              const NodeGeometrySimulationInput *input_data =
+                  static_cast<const NodeGeometrySimulationInput *>(node->storage);
+              LISTBASE_FOREACH (bNode *, output_node, &ntree->nodes) {
+                if (output_node->identifier == input_data->output_node_id) {
+                  /* 'Delta Time' socket in input nodes */
+                  set_socket_identifiers(node, output_node, 1);
+                  break;
+                }
+              }
+            }
+            if (node->type == GEO_NODE_SIMULATION_OUTPUT) {
+              set_socket_identifiers(node, node, 0);
+            }
+          }
+        }
+      }
+    }
   }
 
   /* fcm->name was never used to store modifier name so it has always been an empty string. Now

source/blender/blenloader/intern/versioning_userdef.c

@@ -108,6 +108,7 @@ static void do_versions_theme(const UserDef *userdef, bTheme *btheme)
    */
   {
     /* Keep this block, even when empty. */
+    FROM_DEFAULT_V4_UCHAR(space_node.node_zone_simulation);
   }
 
 #undef FROM_DEFAULT_V4_UCHAR

source/blender/depsgraph/intern/builder/deg_builder_nodes.cc

@@ -87,6 +87,7 @@
 #include "BKE_scene.h"
 #include "BKE_shader_fx.h"
 #include "BKE_simulation.h"
+#include "BKE_simulation_state.hh"
 #include "BKE_sound.h"
 #include "BKE_tracking.h"
 #include "BKE_volume.h"
@@ -906,6 +907,20 @@ void DepsgraphNodeBuilder::build_object_modifiers(Object *object)
   LISTBASE_FOREACH (ModifierData *, modifier, &object->modifiers) {
     OperationNode *modifier_node = add_operation_node(
         &object->id, NodeType::GEOMETRY, OperationCode::MODIFIER, nullptr, modifier->name);
+    if (modifier->type == eModifierType_Nodes) {
+      modifier_node->evaluate = [nmd = reinterpret_cast<NodesModifierData *>(modifier),
+                                 modifier_node](::Depsgraph *depsgraph) {
+        if (!DEG_is_active(depsgraph)) {
+          return;
+        }
+        if (modifier_node->flag & DEPSOP_FLAG_USER_MODIFIED) {
+          if (nmd->simulation_cache &&
+              nmd->simulation_cache->cache_state() == bke::sim::CacheState::Valid) {
+            nmd->simulation_cache->invalidate();
+          }
+        }
+      };
+    }
 
     /* Mute modifier mode if the modifier is not enabled for the dependency graph mode.
      * This handles static (non-animated) mode of the modifier. */

source/blender/editors/include/ED_node.h

@@ -73,6 +73,7 @@ void ED_init_node_socket_type_virtual(struct bNodeSocketType *stype);
 void ED_node_sample_set(const float col[4]);
 void ED_node_draw_snap(
     struct View2D *v2d, const float cent[2], float size, NodeBorder border, unsigned int pos);
+void ED_node_type_draw_color(const char *idname, float *r_color);
 
 /* node_draw.cc */
 

source/blender/editors/include/UI_resources.h

@@ -179,6 +179,8 @@ typedef enum ThemeColorID {
   TH_NODE_GEOMETRY,
   TH_NODE_ATTRIBUTE,
 
+  TH_NODE_ZONE_SIMULATION,
+
   TH_CONSOLE_OUTPUT,
   TH_CONSOLE_INPUT,
   TH_CONSOLE_INFO,

source/blender/editors/interface/resources.cc

@@ -641,6 +641,9 @@ const uchar *UI_ThemeGetColorPtr(bTheme *btheme, int spacetype, int colorid)
         case TH_NODE_GRID_LEVELS:
           cp = &ts->grid_levels;
           break;
+        case TH_NODE_ZONE_SIMULATION:
+          cp = ts->node_zone_simulation;
+          break;
 
         case TH_SEQ_MOVIE:
           cp = ts->movie;

source/blender/editors/object/CMakeLists.txt

@@ -36,6 +36,7 @@ set(SRC
   object_add.cc
   object_bake.cc
   object_bake_api.cc
+  object_bake_simulation.cc
   object_collection.c
   object_constraint.c
   object_data_transfer.c

source/blender/editors/object/object_bake_simulation.cc

@@ -0,0 +1,345 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include <fstream>
+#include <iomanip>
+#include <random>
+
+#include "BLI_endian_defines.h"
+#include "BLI_endian_switch.h"
+#include "BLI_fileops.hh"
+#include "BLI_path_util.h"
+#include "BLI_serialize.hh"
+#include "BLI_vector.hh"
+
+#include "PIL_time.h"
+
+#include "WM_api.h"
+#include "WM_types.h"
+
+#include "ED_screen.h"
+
+#include "DNA_curves_types.h"
+#include "DNA_material_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_pointcloud_types.h"
+#include "DNA_windowmanager_types.h"
+
+#include "BKE_context.h"
+#include "BKE_curves.hh"
+#include "BKE_global.h"
+#include "BKE_instances.hh"
+#include "BKE_lib_id.h"
+#include "BKE_main.h"
+#include "BKE_mesh.hh"
+#include "BKE_object.h"
+#include "BKE_pointcloud.h"
+#include "BKE_report.h"
+#include "BKE_scene.h"
+#include "BKE_simulation_state.hh"
+#include "BKE_simulation_state_serialize.hh"
+
+#include "RNA_access.h"
+#include "RNA_define.h"
+#include "RNA_enum_types.h"
+
+#include "DEG_depsgraph.h"
+#include "DEG_depsgraph_build.h"
+
+#include "object_intern.h"
+
+namespace blender::ed::object::bake_simulation {
+
+static bool bake_simulation_poll(bContext *C)
+{
+  if (!ED_operator_object_active(C)) {
+    return false;
+  }
+  Main *bmain = CTX_data_main(C);
+  const StringRefNull path = BKE_main_blendfile_path(bmain);
+  if (path.is_empty()) {
+    CTX_wm_operator_poll_msg_set(C, "File has to be saved");
+    return false;
+  }
+  return true;
+}
+
+struct ModifierBakeData {
+  NodesModifierData *nmd;
+  std::string meta_dir;
+  std::string bdata_dir;
+  std::unique_ptr<bke::sim::BDataSharing> bdata_sharing;
+};
+
+struct ObjectBakeData {
+  Object *object;
+  Vector<ModifierBakeData> modifiers;
+};
+
+struct BakeSimulationJob {
+  wmWindowManager *wm;
+  Main *bmain;
+  Depsgraph *depsgraph;
+  Scene *scene;
+  Vector<Object *> objects;
+};
+
+static void bake_simulation_job_startjob(void *customdata,
+                                         bool *stop,
+                                         bool *do_update,
+                                         float *progress)
+{
+  using namespace bke::sim;
+
+  BakeSimulationJob &job = *static_cast<BakeSimulationJob *>(customdata);
+  G.is_rendering = true;
+  G.is_break = false;
+  WM_set_locked_interface(job.wm, true);
+
+  Vector<ObjectBakeData> objects_to_bake;
+  for (Object *object : job.objects) {
+    if (!BKE_id_is_editable(job.bmain, &object->id)) {
+      continue;
+    }
+    ObjectBakeData bake_data;
+    bake_data.object = object;
+    LISTBASE_FOREACH (ModifierData *, md, &object->modifiers) {
+      if (md->type == eModifierType_Nodes) {
+        NodesModifierData *nmd = reinterpret_cast<NodesModifierData *>(md);
+        if (nmd->simulation_cache != nullptr) {
+          nmd->simulation_cache->reset();
+        }
+        bake_data.modifiers.append({nmd,
+                                    bke::sim::get_meta_directory(*job.bmain, *object, *md),
+                                    bke::sim::get_bdata_directory(*job.bmain, *object, *md),
+                                    std::make_unique<BDataSharing>()});
+      }
+    }
+    objects_to_bake.append(std::move(bake_data));
+  }
+
+  *progress = 0.0f;
+  *do_update = true;
+
+  const float frame_step_size = 1.0f;
+  const float progress_per_frame = 1.0f / (float(job.scene->r.efra - job.scene->r.sfra + 1) /
+                                           frame_step_size);
+  const int old_frame = job.scene->r.cfra;
+
+  for (float frame_f = job.scene->r.sfra; frame_f <= job.scene->r.efra; frame_f += frame_step_size)
+  {
+    const SubFrame frame{frame_f};
+
+    if (G.is_break || (stop != nullptr && *stop)) {
+      break;
+    }
+
+    job.scene->r.cfra = frame.frame();
+    job.scene->r.subframe = frame.subframe();
+
+    char frame_file_c_str[64];
+    BLI_snprintf(frame_file_c_str, sizeof(frame_file_c_str), "%011.5f", double(frame));
+    BLI_str_replace_char(frame_file_c_str, '.', '_');
+    const StringRefNull frame_file_str = frame_file_c_str;
+
+    BKE_scene_graph_update_for_newframe(job.depsgraph);
+
+    for (ObjectBakeData &object_bake_data : objects_to_bake) {
+      for (ModifierBakeData &modifier_bake_data : object_bake_data.modifiers) {
+        NodesModifierData &nmd = *modifier_bake_data.nmd;
+        if (nmd.simulation_cache == nullptr) {
+          continue;
+        }
+        ModifierSimulationCache &sim_cache = *nmd.simulation_cache;
+        const ModifierSimulationState *sim_state = sim_cache.get_state_at_exact_frame(frame);
+        if (sim_state == nullptr) {
+          continue;
+        }
+
+        const std::string bdata_file_name = frame_file_str + ".bdata";
+        const std::string meta_file_name = frame_file_str + ".json";
+
+        char bdata_path[FILE_MAX];
+        BLI_path_join(bdata_path,
+                      sizeof(bdata_path),
+                      modifier_bake_data.bdata_dir.c_str(),
+                      bdata_file_name.c_str());
+        char meta_path[FILE_MAX];
+        BLI_path_join(meta_path,
+                      sizeof(meta_path),
+                      modifier_bake_data.meta_dir.c_str(),
+                      meta_file_name.c_str());
+
+        BLI_file_ensure_parent_dir_exists(bdata_path);
+        fstream bdata_file{bdata_path, std::ios::out | std::ios::binary};
+        bke::sim::DiskBDataWriter bdata_writer{bdata_file_name, bdata_file, 0};
+
+        io::serialize::DictionaryValue io_root;
+        bke::sim::serialize_modifier_simulation_state(
+            *sim_state, bdata_writer, *modifier_bake_data.bdata_sharing, io_root);
+
+        BLI_file_ensure_parent_dir_exists(meta_path);
+        io::serialize::write_json_file(meta_path, io_root);
+      }
+    }
+
+    *progress += progress_per_frame;
+    *do_update = true;
+  }
+
+  for (ObjectBakeData &object_bake_data : objects_to_bake) {
+    for (ModifierBakeData &modifier_bake_data : object_bake_data.modifiers) {
+      NodesModifierData &nmd = *modifier_bake_data.nmd;
+      if (nmd.simulation_cache) {
+        /* Tag the caches as being baked so that they are not changed anymore. */
+        nmd.simulation_cache->cache_state_ = CacheState::Baked;
+      }
+    }
+    DEG_id_tag_update(&object_bake_data.object->id, ID_RECALC_GEOMETRY);
+  }
+
+  job.scene->r.cfra = old_frame;
+  DEG_time_tag_update(job.bmain);
+
+  *progress = 1.0f;
+  *do_update = true;
+}
+
+static void bake_simulation_job_endjob(void *customdata)
+{
+  BakeSimulationJob &job = *static_cast<BakeSimulationJob *>(customdata);
+  WM_set_locked_interface(job.wm, false);
+  G.is_rendering = false;
+  WM_main_add_notifier(NC_OBJECT | ND_MODIFIER, nullptr);
+}
+
+static int bake_simulation_invoke(bContext *C, wmOperator *op, const wmEvent * /*event*/)
+{
+  wmWindowManager *wm = CTX_wm_manager(C);
+  Scene *scene = CTX_data_scene(C);
+  Depsgraph *depsgraph = CTX_data_depsgraph_pointer(C);
+  Main *bmain = CTX_data_main(C);
+
+  BakeSimulationJob *job = MEM_new<BakeSimulationJob>(__func__);
+  job->wm = wm;
+  job->bmain = bmain;
+  job->depsgraph = depsgraph;
+  job->scene = scene;
+
+  if (RNA_boolean_get(op->ptr, "selected")) {
+    CTX_DATA_BEGIN (C, Object *, object, selected_objects) {
+      job->objects.append(object);
+    }
+    CTX_DATA_END;
+  }
+  else {
+    if (Object *object = CTX_data_active_object(C)) {
+      job->objects.append(object);
+    }
+  }
+
+  wmJob *wm_job = WM_jobs_get(wm,
+                              CTX_wm_window(C),
+                              CTX_data_scene(C),
+                              "Bake Simulation Nodes",
+                              WM_JOB_PROGRESS,
+                              WM_JOB_TYPE_BAKE_SIMULATION_NODES);
+
+  WM_jobs_customdata_set(
+      wm_job, job, [](void *job) { MEM_delete(static_cast<BakeSimulationJob *>(job)); });
+  WM_jobs_timer(wm_job, 0.1, NC_OBJECT | ND_MODIFIER, NC_OBJECT | ND_MODIFIER);
+  WM_jobs_callbacks(
+      wm_job, bake_simulation_job_startjob, nullptr, nullptr, bake_simulation_job_endjob);
+
+  WM_jobs_start(CTX_wm_manager(C), wm_job);
+  WM_event_add_modal_handler(C, op);
+  return OPERATOR_RUNNING_MODAL;
+}
+
+static int bake_simulation_modal(bContext *C, wmOperator * /*op*/, const wmEvent * /*event*/)
+{
+  if (!WM_jobs_test(CTX_wm_manager(C), CTX_data_scene(C), WM_JOB_TYPE_BAKE_SIMULATION_NODES)) {
+    return OPERATOR_FINISHED | OPERATOR_PASS_THROUGH;
+  }
+  return OPERATOR_PASS_THROUGH;
+}
+
+static int delete_baked_simulation_exec(bContext *C, wmOperator *op)
+{
+  Main *bmain = CTX_data_main(C);
+
+  Vector<Object *> objects;
+  if (RNA_boolean_get(op->ptr, "selected")) {
+    CTX_DATA_BEGIN (C, Object *, object, selected_objects) {
+      objects.append(object);
+    }
+    CTX_DATA_END;
+  }
+  else {
+    if (Object *object = CTX_data_active_object(C)) {
+      objects.append(object);
+    }
+  }
+
+  if (objects.is_empty()) {
+    return OPERATOR_CANCELLED;
+  }
+
+  for (Object *object : objects) {
+    LISTBASE_FOREACH (ModifierData *, md, &object->modifiers) {
+      if (md->type == eModifierType_Nodes) {
+        NodesModifierData *nmd = reinterpret_cast<NodesModifierData *>(md);
+        const std::string bake_directory = bke::sim::get_bake_directory(*bmain, *object, *md);
+        if (BLI_exists(bake_directory.c_str())) {
+          if (BLI_delete(bake_directory.c_str(), true, true)) {
+            BKE_reportf(op->reports,
+                        RPT_ERROR,
+                        "Failed to remove bake directory %s",
+                        bake_directory.c_str());
+          }
+        }
+        if (nmd->simulation_cache != nullptr) {
+          nmd->simulation_cache->reset();
+        }
+      }
+    }
+
+    DEG_id_tag_update(&object->id, ID_RECALC_GEOMETRY);
+  }
+
+  WM_event_add_notifier(C, NC_OBJECT | ND_MODIFIER, nullptr);
+
+  return OPERATOR_FINISHED;
+}
+
+}  // namespace blender::ed::object::bake_simulation
+
+void OBJECT_OT_simulation_nodes_cache_bake(wmOperatorType *ot)
+{
+  using namespace blender::ed::object::bake_simulation;
+
+  ot->name = "Bake Simulation";
+  ot->description = "Bake simulations in geometry nodes modifiers";
+  ot->idname = __func__;
+
+  ot->invoke = bake_simulation_invoke;
+  ot->modal = bake_simulation_modal;
+  ot->poll = bake_simulation_poll;
+
+  RNA_def_boolean(ot->srna, "selected", false, "Selected", "Bake cache on all selected objects");
+}
+
+void OBJECT_OT_simulation_nodes_cache_delete(wmOperatorType *ot)
+{
+  using namespace blender::ed::object::bake_simulation;
+
+  ot->name = "Delete Cached Simulation";
+  ot->description = "Delete cached/baked simulations in geometry nodes modifiers";
+  ot->idname = __func__;
+
+  ot->exec = delete_baked_simulation_exec;
+  ot->poll = ED_operator_object_active;
+
+  RNA_def_boolean(ot->srna, "selected", false, "Selected", "Delete cache on all selected objects");
+}

source/blender/editors/object/object_intern.h

@@ -338,6 +338,11 @@ void OBJECT_OT_collection_objects_select(struct wmOperatorType *ot);
 void OBJECT_OT_bake_image(wmOperatorType *ot);
 void OBJECT_OT_bake(wmOperatorType *ot);
 
+/* object_bake_simulation.cc */
+
+void OBJECT_OT_simulation_nodes_cache_bake(wmOperatorType *ot);
+void OBJECT_OT_simulation_nodes_cache_delete(wmOperatorType *ot);
+
 /* object_random.c */
 
 void TRANSFORM_OT_vertex_random(struct wmOperatorType *ot);

source/blender/editors/object/object_ops.c

@@ -260,6 +260,8 @@ void ED_operatortypes_object(void)
 
   WM_operatortype_append(OBJECT_OT_bake_image);
   WM_operatortype_append(OBJECT_OT_bake);
+  WM_operatortype_append(OBJECT_OT_simulation_nodes_cache_bake);
+  WM_operatortype_append(OBJECT_OT_simulation_nodes_cache_delete);
   WM_operatortype_append(OBJECT_OT_drop_named_material);
   WM_operatortype_append(OBJECT_OT_drop_geometry_nodes);
   WM_operatortype_append(OBJECT_OT_unlink_data);

source/blender/editors/space_action/action_draw.cc

@@ -21,6 +21,7 @@
 #include "DNA_anim_types.h"
 #include "DNA_cachefile_types.h"
 #include "DNA_gpencil_legacy_types.h"
+#include "DNA_modifier_types.h"
 #include "DNA_object_types.h"
 #include "DNA_scene_types.h"
 #include "DNA_screen_types.h"
@@ -28,6 +29,7 @@
 #include "BKE_action.h"
 #include "BKE_context.h"
 #include "BKE_pointcache.h"
+#include "BKE_simulation_state.hh"
 
 /* Everything from source (BIF, BDR, BSE) ------------------------------ */
 
@@ -679,6 +681,51 @@ static void timeline_cache_draw_single(PTCacheID *pid, float y_offset, float hei
   GPU_matrix_pop();
 }
 
+static void timeline_cache_draw_simulation_nodes(
+    const Scene &scene,
+    const blender::bke::sim::ModifierSimulationCache &cache,
+    const float y_offset,
+    const float height,
+    const uint pos_id)
+{
+  GPU_matrix_push();
+  GPU_matrix_translate_2f(0.0, (float)V2D_SCROLL_HANDLE_HEIGHT + y_offset);
+  GPU_matrix_scale_2f(1.0, height);
+
+  float color[4];
+  switch (cache.cache_state()) {
+    case blender::bke::sim::CacheState::Invalid: {
+      copy_v4_fl4(color, 0.8, 0.8, 0.2, 0.3);
+      break;
+    }
+    case blender::bke::sim::CacheState::Valid: {
+      copy_v4_fl4(color, 0.8, 0.8, 0.2, 1.0);
+      break;
+    }
+    case blender::bke::sim::CacheState::Baked: {
+      copy_v4_fl4(color, 1.0, 0.6, 0.2, 1.0);
+      break;
+    }
+  }
+
+  immUniformColor4fv(color);
+
+  const int start_frame = scene.r.sfra;
+  const int end_frame = scene.r.efra;
+  const int frames_num = end_frame - start_frame + 1;
+  const blender::IndexRange frames_range(start_frame, frames_num);
+
+  immBeginAtMost(GPU_PRIM_TRIS, frames_num * 6);
+  for (const int frame : frames_range) {
+    if (cache.has_state_at_frame(frame)) {
+      immRectf_fast(pos_id, frame - 0.5f, 0, frame + 0.5f, 1.0f);
+    }
+  }
+  immEnd();
+
+  GPU_matrix_pop();
+}
+
 void timeline_draw_cache(const SpaceAction *saction, const Object *ob, const Scene *scene)
 {
   if ((saction->cache_display & TIME_CACHE_DISPLAY) == 0 || ob == nullptr) {
@@ -710,6 +757,16 @@ void timeline_draw_cache(const SpaceAction *saction, const Object *ob, const Sce
 
     y_offset += cache_draw_height;
   }
+  LISTBASE_FOREACH (ModifierData *, md, &ob->modifiers) {
+    if (md->type == eModifierType_Nodes) {
+      const NodesModifierData *nmd = reinterpret_cast<NodesModifierData *>(md);
+      if (nmd->simulation_cache != nullptr) {
+        timeline_cache_draw_simulation_nodes(
+            *scene, *nmd->simulation_cache, y_offset, cache_draw_height, pos_id);
+        y_offset += cache_draw_height;
+      }
+    }
+  }
 
   GPU_blend(GPU_BLEND_NONE);
   immUnbindProgram();

source/blender/editors/space_node/CMakeLists.txt

@@ -12,6 +12,7 @@ set(INC
   ../../depsgraph
   ../../draw
   ../../functions
+  ../../geometry
   ../../gpu
   ../../imbuf
   ../../makesdna

source/blender/editors/space_node/clipboard.cc

@@ -16,6 +16,8 @@
 #include "ED_render.h"
 #include "ED_screen.h"
 
+#include "NOD_socket.h"
+
 #include "RNA_access.h"
 #include "RNA_define.h"
 
@@ -182,6 +184,33 @@ void NODE_OT_clipboard_copy(wmOperatorType *ot)
 /** \name Paste
  * \{ */
 
+static void remap_pairing(bNodeTree &dst_tree, const Map<const bNode *, bNode *> &node_map)
+{
+  /* We don't have the old tree for looking up output nodes by ID,
+   * so we have to build a map first to find copied output nodes in the new tree. */
+  Map<int32_t, bNode *> dst_output_node_map;
+  for (const auto &item : node_map.items()) {
+    if (item.key->type == GEO_NODE_SIMULATION_OUTPUT) {
+      dst_output_node_map.add_new(item.key->identifier, item.value);
+    }
+  }
+
+  for (bNode *dst_node : node_map.values()) {
+    if (dst_node->type == GEO_NODE_SIMULATION_INPUT) {
+      NodeGeometrySimulationInput &data = *static_cast<NodeGeometrySimulationInput *>(
+          dst_node->storage);
+      if (const bNode *output_node = dst_output_node_map.lookup_default(data.output_node_id,
+                                                                        nullptr)) {
+        data.output_node_id = output_node->identifier;
+      }
+      else {
+        data.output_node_id = 0;
+        blender::nodes::update_node_declaration_and_sockets(dst_tree, *dst_node);
+      }
+    }
+  }
+}
+
 static int node_clipboard_paste_exec(bContext *C, wmOperator *op)
 {
   SpaceNode &snode = *CTX_wm_space_node(C);
@@ -285,6 +314,12 @@ static int node_clipboard_paste_exec(bContext *C, wmOperator *op)
     }
   }
 
+  for (bNode *new_node : node_map.values()) {
+    nodeDeclarationEnsure(&tree, new_node);
+  }
+
+  remap_pairing(tree, node_map);
+
   tree.ensure_topology_cache();
   for (bNode *new_node : node_map.values()) {
     /* Update multi input socket indices in case all connected nodes weren't copied. */

source/blender/editors/space_node/drawnode.cc

@@ -1498,6 +1498,23 @@ void ED_init_node_socket_type_virtual(bNodeSocketType *stype)
   stype->draw_color = node_socket_virtual_draw_color;
 }
 
+void ED_node_type_draw_color(const char *idname, float *r_color)
+{
+  using namespace blender::ed::space_node;
+
+  const bNodeSocketType *typeinfo = nodeSocketTypeFind(idname);
+  if (!typeinfo || typeinfo->type == SOCK_CUSTOM) {
+    r_color[0] = 0.0f;
+    r_color[1] = 0.0f;
+    r_color[2] = 0.0f;
+    r_color[3] = 0.0f;
+    return;
+  }
+
+  BLI_assert(typeinfo->type < ARRAY_SIZE(std_node_socket_colors));
+  copy_v4_v4(r_color, std_node_socket_colors[typeinfo->type]);
+}
+
 namespace blender::ed::space_node {
 
 /* ************** Generic drawing ************** */
@@ -2220,8 +2237,8 @@ void node_draw_link(const bContext &C,
   node_draw_link_bezier(C, v2d, snode, link, th_col1, th_col2, th_col3, selected);
 }
 
-static std::array<float2, 4> node_link_bezier_points_dragged(const SpaceNode &snode,
-                                                             const bNodeLink &link)
+std::array<float2, 4> node_link_bezier_points_dragged(const SpaceNode &snode,
+                                                      const bNodeLink &link)
 {
   const float2 cursor = snode.runtime->cursor * UI_SCALE_FAC;
   std::array<float2, 4> points;

source/blender/editors/space_node/node_draw.cc

@@ -22,6 +22,7 @@
 #include "DNA_world_types.h"
 
 #include "BLI_array.hh"
+#include "BLI_convexhull_2d.h"
 #include "BLI_map.hh"
 #include "BLI_set.hh"
 #include "BLI_span.hh"
@@ -32,12 +33,15 @@
 
 #include "BKE_compute_contexts.hh"
 #include "BKE_context.h"
+#include "BKE_curves.hh"
+#include "BKE_global.h"
 #include "BKE_idtype.h"
 #include "BKE_lib_id.h"
 #include "BKE_main.h"
 #include "BKE_node.h"
 #include "BKE_node_runtime.hh"
 #include "BKE_node_tree_update.h"
+#include "BKE_node_tree_zones.hh"
 #include "BKE_object.h"
 #include "BKE_type_conversions.hh"
 
@@ -80,11 +84,15 @@
 #include "FN_field.hh"
 #include "FN_field_cpp_type.hh"
 
+#include "GEO_fillet_curves.hh"
+
 #include "../interface/interface_intern.hh" /* TODO: Remove */
 
 #include "node_intern.hh" /* own include */
 
 namespace geo_log = blender::nodes::geo_eval_log;
+using blender::bke::node_tree_zones::TreeZone;
+using blender::bke::node_tree_zones::TreeZones;
 
 /**
  * This is passed to many functions which draw the node editor.
@@ -2149,7 +2157,9 @@ static void node_draw_basis(const bContext &C,
   }
 
   /* Shadow. */
-  node_draw_shadow(snode, node, BASIS_RAD, 1.0f);
+  if (!ELEM(node.type, GEO_NODE_SIMULATION_INPUT, GEO_NODE_SIMULATION_OUTPUT)) {
+    node_draw_shadow(snode, node, BASIS_RAD, 1.0f);
+  }
 
   const rctf &rct = node.runtime->totr;
   float color[4];
@@ -2423,6 +2433,10 @@ static void node_draw_basis(const bContext &C,
     else if (nodeTypeUndefined(&node)) {
       UI_GetThemeColor4fv(TH_REDALERT, color_outline);
     }
+    else if (ELEM(node.type, GEO_NODE_SIMULATION_INPUT, GEO_NODE_SIMULATION_OUTPUT)) {
+      UI_GetThemeColor4fv(TH_NODE_ZONE_SIMULATION, color_outline);
+      color_outline[3] = 1.0f;
+    }
     else {
       UI_GetThemeColorBlendShade4fv(TH_BACK, TH_NODE, 0.4f, -20, color_outline);
     }
@@ -3036,6 +3050,173 @@ static void node_draw(const bContext &C,
   }
 }
 
+static void add_rect_corner_positions(Vector<float2> &positions, const rctf &rect)
+{
+  positions.append({rect.xmin, rect.ymin});
+  positions.append({rect.xmin, rect.ymax});
+  positions.append({rect.xmax, rect.ymin});
+  positions.append({rect.xmax, rect.ymax});
+}
+
+static void find_bounds_by_zone_recursive(const SpaceNode &snode,
+                                          const TreeZone &zone,
+                                          const Span<std::unique_ptr<TreeZone>> all_zones,
+                                          MutableSpan<Vector<float2>> r_bounds_by_zone)
+{
+  const float node_padding = UI_UNIT_X;
+  const float zone_padding = 0.3f * UI_UNIT_X;
+
+  Vector<float2> &bounds = r_bounds_by_zone[zone.index];
+  if (!bounds.is_empty()) {
+    return;
+  }
+
+  Vector<float2> possible_bounds;
+  for (const TreeZone *child_zone : zone.child_zones) {
+    find_bounds_by_zone_recursive(snode, *child_zone, all_zones, r_bounds_by_zone);
+    const Span<float2> child_bounds = r_bounds_by_zone[child_zone->index];
+    for (const float2 &pos : child_bounds) {
+      rctf rect;
+      BLI_rctf_init_pt_radius(&rect, pos, zone_padding);
+      add_rect_corner_positions(possible_bounds, rect);
+    }
+  }
+  for (const bNode *child_node : zone.child_nodes) {
+    rctf rect = child_node->runtime->totr;
+    BLI_rctf_pad(&rect, node_padding, node_padding);
+    add_rect_corner_positions(possible_bounds, rect);
+  }
+  if (zone.input_node) {
+    const rctf &totr = zone.input_node->runtime->totr;
+    rctf rect = totr;
+    BLI_rctf_pad(&rect, node_padding, node_padding);
+    rect.xmin = math::interpolate(totr.xmin, totr.xmax, 0.25f);
+    add_rect_corner_positions(possible_bounds, rect);
+  }
+  if (zone.output_node) {
+    const rctf &totr = zone.output_node->runtime->totr;
+    rctf rect = totr;
+    BLI_rctf_pad(&rect, node_padding, node_padding);
+    rect.xmax = math::interpolate(totr.xmin, totr.xmax, 0.75f);
+    add_rect_corner_positions(possible_bounds, rect);
+  }
+
+  if (snode.runtime->linkdrag) {
+    for (const bNodeLink &link : snode.runtime->linkdrag->links) {
+      if (link.fromnode == nullptr) {
+        continue;
+      }
+      if (zone.contains_node_recursively(*link.fromnode) || zone.input_node == link.fromnode) {
+        const float2 pos = node_link_bezier_points_dragged(snode, link)[3];
+        rctf rect;
+        BLI_rctf_init_pt_radius(&rect, pos, node_padding);
+        add_rect_corner_positions(possible_bounds, rect);
+      }
+    }
+  }
+
+  Vector<int> convex_indices(possible_bounds.size());
+  const int convex_positions_num = BLI_convexhull_2d(
+      reinterpret_cast<float(*)[2]>(possible_bounds.data()),
+      possible_bounds.size(),
+      convex_indices.data());
+  convex_indices.resize(convex_positions_num);
+
+  for (const int i : convex_indices) {
+    bounds.append(possible_bounds[i]);
+  }
+}
+
+static void node_draw_zones(TreeDrawContext & /*tree_draw_ctx*/,
+                            const ARegion &region,
+                            const SpaceNode &snode,
+                            const bNodeTree &ntree)
+{
+  const TreeZones *zones = bke::node_tree_zones::get_tree_zones(ntree);
+  if (zones == nullptr) {
+    return;
+  }
+
+  Array<Vector<float2>> bounds_by_zone(zones->zones.size());
+  Array<bke::CurvesGeometry> fillet_curve_by_zone(zones->zones.size());
+
+  for (const int zone_i : zones->zones.index_range()) {
+    const TreeZone &zone = *zones->zones[zone_i];
+
+    find_bounds_by_zone_recursive(snode, zone, zones->zones, bounds_by_zone);
+    const Span<float2> boundary_positions = bounds_by_zone[zone_i];
+    const int boundary_positions_num = boundary_positions.size();
+
+    bke::CurvesGeometry boundary_curve(boundary_positions_num, 1);
+    boundary_curve.cyclic_for_write().first() = true;
+    boundary_curve.fill_curve_types(CURVE_TYPE_POLY);
+    MutableSpan<float3> boundary_curve_positions = boundary_curve.positions_for_write();
+    boundary_curve.offsets_for_write().copy_from({0, boundary_positions_num});
+    for (const int i : boundary_positions.index_range()) {
+      boundary_curve_positions[i] = float3(boundary_positions[i], 0.0f);
+    }
+
+    fillet_curve_by_zone[zone_i] = geometry::fillet_curves_poly(
+        boundary_curve,
+        IndexRange(1),
+        VArray<float>::ForSingle(BASIS_RAD, boundary_positions_num),
+        VArray<int>::ForSingle(5, boundary_positions_num),
+        true,
+        {});
+  }
+
+  const View2D &v2d = region.v2d;
+  float scale;
+  UI_view2d_scale_get(&v2d, &scale, nullptr);
+  float line_width = 1.0f * scale;
+  float viewport[4] = {};
+  GPU_viewport_size_get_f(viewport);
+  float zone_color[4];
+  UI_GetThemeColor4fv(TH_NODE_ZONE_SIMULATION, zone_color);
+
+  const uint pos = GPU_vertformat_attr_add(
+      immVertexFormat(), "pos", GPU_COMP_F32, 3, GPU_FETCH_FLOAT);
+
+  /* Draw all the contour lines after to prevent them from getting hidden by overlapping zones. */
+  for (const int zone_i : zones->zones.index_range()) {
+    if (zone_color[3] == 0.0f) {
+      break;
+    }
+    const Span<float3> fillet_boundary_positions = fillet_curve_by_zone[zone_i].positions();
+    /* Draw the background. */
+    immBindBuiltinProgram(GPU_SHADER_3D_UNIFORM_COLOR);
+    immUniformThemeColorBlend(TH_BACK, TH_NODE_ZONE_SIMULATION, zone_color[3]);
+
+    immBegin(GPU_PRIM_TRI_FAN, fillet_boundary_positions.size() + 1);
+    for (const float3 &p : fillet_boundary_positions) {
+      immVertex3fv(pos, p);
+    }
+    immVertex3fv(pos, fillet_boundary_positions[0]);
+    immEnd();
+
+    immUnbindProgram();
+  }
+
+  for (const int zone_i : zones->zones.index_range()) {
+    const Span<float3> fillet_boundary_positions = fillet_curve_by_zone[zone_i].positions();
+    /* Draw the contour lines. */
+    immBindBuiltinProgram(GPU_SHADER_3D_POLYLINE_UNIFORM_COLOR);
+
+    immUniform2fv("viewportSize", &viewport[2]);
+    immUniform1f("lineWidth", line_width * U.pixelsize);
+
+    immUniformThemeColorAlpha(TH_NODE_ZONE_SIMULATION, 1.0f);
+    immBegin(GPU_PRIM_LINE_STRIP, fillet_boundary_positions.size() + 1);
+    for (const float3 &p : fillet_boundary_positions) {
+      immVertex3fv(pos, p);
+    }
+    immVertex3fv(pos, fillet_boundary_positions[0]);
+    immEnd();
+
+    immUnbindProgram();
+  }
+}
+
 #define USE_DRAW_TOT_UPDATE
 
 static void node_draw_nodetree(const bContext &C,
@@ -3203,6 +3384,7 @@ static void draw_nodetree(const bContext &C,
   }
 
   node_update_nodetree(C, tree_draw_ctx, ntree, nodes, blocks);
+  node_draw_zones(tree_draw_ctx, region, *snode, ntree);
   node_draw_nodetree(C, tree_draw_ctx, region, *snode, ntree, nodes, blocks, parent_key);
 }
 
@@ -3229,6 +3411,9 @@ static void draw_background_color(const SpaceNode &snode)
 
 void node_draw_space(const bContext &C, ARegion &region)
 {
+  if (G.is_rendering) {
+    return;
+  }
   wmWindow *win = CTX_wm_window(&C);
   SpaceNode &snode = *CTX_wm_space_node(&C);
   View2D &v2d = region.v2d;

source/blender/editors/space_node/node_edit.cc

@@ -65,6 +65,7 @@
 #include "NOD_composite.h"
 #include "NOD_geometry.h"
 #include "NOD_shader.h"
+#include "NOD_socket.h"
 #include "NOD_texture.h"
 #include "node_intern.hh" /* own include */
 
@@ -1250,6 +1251,33 @@ static void node_duplicate_reparent_recursive(bNodeTree *ntree,
   }
 }
 
+static void remap_pairing(bNodeTree &dst_tree, const Map<bNode *, bNode *> &node_map)
+{
+  /* We don't have the old tree for looking up output nodes by ID,
+   * so we have to build a map first to find copied output nodes in the new tree. */
+  Map<uint32_t, bNode *> dst_output_node_map;
+  for (const auto &item : node_map.items()) {
+    if (item.key->type == GEO_NODE_SIMULATION_OUTPUT) {
+      dst_output_node_map.add_new(item.key->identifier, item.value);
+    }
+  }
+
+  for (bNode *dst_node : node_map.values()) {
+    if (dst_node->type == GEO_NODE_SIMULATION_INPUT) {
+      NodeGeometrySimulationInput *data = static_cast<NodeGeometrySimulationInput *>(
+          dst_node->storage);
+      if (const bNode *output_node = dst_output_node_map.lookup_default(data->output_node_id,
+                                                                        nullptr)) {
+        data->output_node_id = output_node->identifier;
+      }
+      else {
+        data->output_node_id = 0;
+        blender::nodes::update_node_declaration_and_sockets(dst_tree, *dst_node);
+      }
+    }
+  }
+}
+
 static int node_duplicate_exec(bContext *C, wmOperator *op)
 {
   Main *bmain = CTX_data_main(C);
@@ -1323,6 +1351,10 @@ static int node_duplicate_exec(bContext *C, wmOperator *op)
     }
   }
 
+  for (bNode *node : node_map.values()) {
+    nodeDeclarationEnsure(ntree, node);
+  }
+
   /* Clear flags for recursive depth-first iteration. */
   for (bNode *node : ntree->all_nodes()) {
     node->flag &= ~NODE_TEST;
@@ -1334,6 +1366,8 @@ static int node_duplicate_exec(bContext *C, wmOperator *op)
     }
   }
 
+  remap_pairing(*ntree, node_map);
+
   /* Deselect old nodes, select the copies instead. */
   for (const auto item : node_map.items()) {
     bNode *src_node = item.key;
@@ -1768,6 +1802,9 @@ static int node_delete_exec(bContext *C, wmOperator * /*op*/)
 
   ED_preview_kill_jobs(CTX_wm_manager(C), bmain);
 
+  /* Delete paired nodes as well. */
+  node_select_paired(*snode->edittree);
+
   LISTBASE_FOREACH_MUTABLE (bNode *, node, &snode->edittree->nodes) {
     if (node->flag & SELECT) {
       nodeRemoveNode(bmain, snode->edittree, node, true);
@@ -1855,6 +1892,9 @@ static int node_delete_reconnect_exec(bContext *C, wmOperator * /*op*/)
 
   ED_preview_kill_jobs(CTX_wm_manager(C), CTX_data_main(C));
 
+  /* Delete paired nodes as well. */
+  node_select_paired(*snode->edittree);
+
   LISTBASE_FOREACH_MUTABLE (bNode *, node, &snode->edittree->nodes) {
     if (node->flag & SELECT) {
       nodeInternalRelink(snode->edittree, node);

source/blender/editors/space_node/node_group.cc

@@ -440,6 +440,28 @@ void NODE_OT_group_ungroup(wmOperatorType *ot)
 /** \name Separate Operator
  * \{ */
 
+static void remap_pairing(bNodeTree &dst_tree, const Set<bNode *> &nodes)
+{
+  for (bNode *dst_node : nodes) {
+    if (dst_node->type == GEO_NODE_SIMULATION_INPUT) {
+      NodeGeometrySimulationInput &data = *static_cast<NodeGeometrySimulationInput *>(
+          dst_node->storage);
+      /* XXX Technically this is not correct because the output_node_id is only valid
+       * in the original node group tree and we'd have map old IDs to new nodes first.
+       * The ungroup operator does not build a node map, it just expects node IDs to
+       * remain unchanged, which is probably true most of the time because they are
+       * mostly random numbers out of the uint32_t range. */
+      if (const bNode *output_node = dst_tree.node_by_id(data.output_node_id)) {
+        data.output_node_id = output_node->identifier;
+      }
+      else {
+        data.output_node_id = 0;
+        blender::nodes::update_node_declaration_and_sockets(dst_tree, *dst_node);
+      }
+    }
+  }
+}
+
 /**
  * \return True if successful.
  */
@@ -457,10 +479,13 @@ static bool node_group_separate_selected(
   nodes_to_move.remove_if(
       [](const bNode *node) { return node->is_group_input() || node->is_group_output(); });
 
+  Set<bNode *> new_nodes;
+
   for (bNode *node : nodes_to_move) {
     bNode *newnode;
     if (make_copy) {
       newnode = bke::node_copy_with_mapping(&ntree, *node, LIB_ID_COPY_DEFAULT, true, socket_map);
+      new_nodes.add_new(newnode);
     }
     else {
       newnode = node;
@@ -525,6 +550,10 @@ static bool node_group_separate_selected(
     }
   }
 
+  for (bNode *node : new_nodes) {
+    nodeDeclarationEnsure(&ntree, node);
+  }
+
   /* and copy across the animation,
    * note that the animation data's action can be nullptr here */
   if (ngroup.adt) {
@@ -537,6 +566,8 @@ static bool node_group_separate_selected(
     }
   }
 
+  remap_pairing(ntree, new_nodes);
+
   BKE_ntree_update_tag_all(&ntree);
   if (!make_copy) {
     BKE_ntree_update_tag_all(&ngroup);

source/blender/editors/space_node/node_intern.hh

@@ -190,6 +190,10 @@ void node_socket_select(bNode *node, bNodeSocket &sock);
 void node_socket_deselect(bNode *node, bNodeSocket &sock, bool deselect_node);
 void node_deselect_all_input_sockets(bNodeTree &node_tree, bool deselect_nodes);
 void node_deselect_all_output_sockets(bNodeTree &node_tree, bool deselect_nodes);
+/**
+ * Select nodes that are paired to a selected node.
+ */
+void node_select_paired(bNodeTree &node_tree);
 void node_select_single(bContext &C, bNode &node);
 
 void NODE_OT_select(wmOperatorType *ot);
@@ -247,6 +251,8 @@ void node_draw_link_bezier(const bContext &C,
                            int th_col3,
                            bool selected);
 
+std::array<float2, 4> node_link_bezier_points_dragged(const SpaceNode &snode,
+                                                      const bNodeLink &link);
 void node_link_bezier_points_evaluated(const bNodeLink &link,
                                        std::array<float2, NODE_LINK_RESOL + 1> &coords);
 

source/blender/editors/space_node/node_select.cc

@@ -309,6 +309,21 @@ void node_deselect_all_output_sockets(bNodeTree &node_tree, const bool deselect_
   }
 }
 
+void node_select_paired(bNodeTree &node_tree)
+{
+  for (bNode *input_node : node_tree.nodes_by_type("GeometryNodeSimulationInput")) {
+    const auto *storage = static_cast<const NodeGeometrySimulationInput *>(input_node->storage);
+    if (bNode *output_node = node_tree.node_by_id(storage->output_node_id)) {
+      if (input_node->flag & NODE_SELECT) {
+        output_node->flag |= NODE_SELECT;
+      }
+      if (output_node->flag & NODE_SELECT) {
+        input_node->flag |= NODE_SELECT;
+      }
+    }
+  }
+}
+
 VectorSet<bNode *> get_selected_nodes(bNodeTree &node_tree)
 {
   VectorSet<bNode *> selected_nodes;

source/blender/makesdna/DNA_modifier_types.h

@@ -10,6 +10,15 @@
 #include "DNA_listBase.h"
 #include "DNA_session_uuid_types.h"
 
+#ifdef __cplusplus
+namespace blender::bke::sim {
+class ModifierSimulationCache;
+}
+using ModifierSimulationCacheHandle = blender::bke::sim::ModifierSimulationCache;
+#else
+typedef struct ModifierSimulationCacheHandle ModifierSimulationCacheHandle;
+#endif
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -2316,7 +2325,8 @@ typedef struct NodesModifierData {
    * This can be used to help the user to debug a node tree.
    */
   void *runtime_eval_log;
-  void *_pad1;
+
+  ModifierSimulationCacheHandle *simulation_cache;
 } NodesModifierData;
 
 typedef struct MeshToVolumeModifierData {

source/blender/makesdna/DNA_node_types.h

@@ -16,6 +16,8 @@
 #ifdef __cplusplus
 namespace blender {
 template<typename T> class Span;
+template<typename T> class MutableSpan;
+class IndexRange;
 class StringRef;
 class StringRefNull;
 }  // namespace blender
@@ -1594,6 +1596,39 @@ typedef struct NodeGeometryUVUnwrap {
   uint8_t method;
 } NodeGeometryUVUnwrap;
 
+typedef struct NodeSimulationItem {
+  char *name;
+  /** #eNodeSocketDatatype. */
+  short socket_type;
+  /** #eAttrDomain. */
+  short attribute_domain;
+  /**
+   * Generates unique identifier for sockets which stays the same even when the item order or
+   * names change.
+   */
+  int identifier;
+} NodeSimulationItem;
+
+typedef struct NodeGeometrySimulationInput {
+  /** bNode.identifier of the corresponding output node. */
+  int32_t output_node_id;
+} NodeGeometrySimulationInput;
+
+typedef struct NodeGeometrySimulationOutput {
+  NodeSimulationItem *items;
+  int items_num;
+  int active_index;
+  /** Number to give unique IDs to state items. */
+  int next_identifier;
+  int _pad;
+
+#ifdef __cplusplus
+  blender::Span<NodeSimulationItem> items_span() const;
+  blender::MutableSpan<NodeSimulationItem> items_span_for_write();
+  blender::IndexRange items_range() const;
+#endif
+} NodeGeometrySimulationOutput;
+
 typedef struct NodeGeometryDistributePointsInVolume {
   /* GeometryNodePointDistributeVolumeMode. */
   uint8_t mode;

source/blender/makesdna/DNA_userdef_types.h

@@ -344,6 +344,9 @@ typedef struct ThemeSpace {
   unsigned char nodeclass_pattern[4], nodeclass_layout[4];
   unsigned char nodeclass_geometry[4], nodeclass_attribute[4];
 
+  unsigned char node_zone_simulation[4];
+  char _pad8[4];
+
   /** For sequence editor. */
   unsigned char movie[4], movieclip[4], mask[4], image[4], scene[4], audio[4];
   unsigned char effect[4], transition[4], meta[4], text_strip[4], color_strip[4];

source/blender/makesrna/RNA_access.h

@@ -216,6 +216,9 @@ bool RNA_enum_name(const EnumPropertyItem *item, int value, const char **r_name)
 bool RNA_enum_description(const EnumPropertyItem *item, int value, const char **description);
 int RNA_enum_from_value(const EnumPropertyItem *item, int value);
 int RNA_enum_from_identifier(const EnumPropertyItem *item, const char *identifier);
+bool RNA_enum_value_from_identifier(const EnumPropertyItem *item,
+                                    const char *identifier,
+                                    int *r_value);
 /**
  * Take care using this with translated enums,
  * prefer #RNA_enum_from_identifier where possible.

source/blender/makesrna/intern/rna_access.cc

@@ -1761,6 +1761,18 @@ int RNA_enum_from_identifier(const EnumPropertyItem *item, const char *identifie
   return -1;
 }
 
+bool RNA_enum_value_from_identifier(const EnumPropertyItem *item,
+                                    const char *identifier,
+                                    int *r_value)
+{
+  const int i = RNA_enum_from_identifier(item, identifier);
+  if (i == -1) {
+    return false;
+  }
+  *r_value = item[i].value;
+  return true;
+}
+
 int RNA_enum_from_name(const EnumPropertyItem *item, const char *name)
 {
   int i = 0;

source/blender/makesrna/intern/rna_nodetree.c

@@ -52,6 +52,8 @@
 #include "RE_texture.h"
 
 #include "NOD_composite.h"
+#include "NOD_geometry.h"
+#include "NOD_socket.h"
 
 #include "DEG_depsgraph.h"
 #include "DEG_depsgraph_query.h"
@@ -4090,6 +4092,167 @@ static void rna_NodeCryptomatte_update_remove(Main *bmain, Scene *scene, Pointer
   rna_Node_update(bmain, scene, ptr);
 }
 
+static void rna_SimulationStateItem_update(Main *bmain, Scene *UNUSED(scene), PointerRNA *ptr)
+{
+  bNodeTree *ntree = (bNodeTree *)ptr->owner_id;
+  NodeSimulationItem *item = (NodeSimulationItem *)ptr->data;
+  bNode *node = NOD_geometry_simulation_output_find_node_by_item(ntree, item);
+
+  BKE_ntree_update_tag_node_property(ntree, node);
+  ED_node_tree_propagate_change(NULL, bmain, ntree);
+}
+
+static bool rna_SimulationStateItem_socket_type_supported(const EnumPropertyItem *item)
+{
+  return NOD_geometry_simulation_output_item_socket_type_supported(
+      (eNodeSocketDatatype)item->value);
+}
+
+static const EnumPropertyItem *rna_SimulationStateItem_socket_type_itemf(bContext *UNUSED(C),
+                                                                         PointerRNA *UNUSED(ptr),
+                                                                         PropertyRNA *UNUSED(prop),
+                                                                         bool *r_free)
+{
+  *r_free = true;
+  return itemf_function_check(node_socket_data_type_items,
+                              rna_SimulationStateItem_socket_type_supported);
+}
+
+static void rna_SimulationStateItem_name_set(PointerRNA *ptr, const char *value)
+{
+  bNodeTree *ntree = (bNodeTree *)ptr->owner_id;
+  NodeSimulationItem *item = (NodeSimulationItem *)ptr->data;
+  bNode *node = NOD_geometry_simulation_output_find_node_by_item(ntree, item);
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+
+  const char *defname = nodeStaticSocketLabel(item->socket_type, 0);
+  NOD_geometry_simulation_output_item_set_unique_name(sim, item, value, defname);
+}
+
+static void rna_SimulationStateItem_color_get(PointerRNA *ptr, float *values)
+{
+  NodeSimulationItem *item = (NodeSimulationItem *)ptr->data;
+
+  const char *socket_type_idname = nodeStaticSocketType(item->socket_type, 0);
+  ED_node_type_draw_color(socket_type_idname, values);
+}
+
+static PointerRNA rna_NodeGeometrySimulationInput_paired_output_get(PointerRNA *ptr)
+{
+  bNodeTree *ntree = (bNodeTree *)ptr->owner_id;
+  bNode *node = (bNode *)ptr->data;
+  bNode *output_node = NOD_geometry_simulation_input_get_paired_output(ntree, node);
+  PointerRNA r_ptr;
+  RNA_pointer_create(&ntree->id, &RNA_Node, output_node, &r_ptr);
+  return r_ptr;
+}
+
+static bool rna_GeometryNodeSimulationInput_pair_with_output(
+    ID *id, bNode *node, bContext *C, ReportList *reports, bNode *output_node)
+{
+  bNodeTree *ntree = (bNodeTree *)id;
+
+  if (!NOD_geometry_simulation_input_pair_with_output(ntree, node, output_node)) {
+    BKE_reportf(reports,
+                RPT_ERROR,
+                "Failed to pair simulation input node %s with output node %s",
+                node->name,
+                output_node->name);
+    return false;
+  }
+
+  BKE_ntree_update_tag_node_property(ntree, node);
+  ED_node_tree_propagate_change(C, CTX_data_main(C), ntree);
+  WM_main_add_notifier(NC_NODE | NA_EDITED, ntree);
+
+  return true;
+}
+
+static NodeSimulationItem *rna_NodeGeometrySimulationOutput_items_new(
+    ID *id, bNode *node, Main *bmain, ReportList *reports, int socket_type, const char *name)
+{
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+  NodeSimulationItem *item = NOD_geometry_simulation_output_add_item(
+      sim, (short)socket_type, name);
+
+  if (item == NULL) {
+    BKE_report(reports, RPT_ERROR, "Unable to create socket");
+  }
+  else {
+    bNodeTree *ntree = (bNodeTree *)id;
+    BKE_ntree_update_tag_node_property(ntree, node);
+    ED_node_tree_propagate_change(NULL, bmain, ntree);
+    WM_main_add_notifier(NC_NODE | NA_EDITED, ntree);
+  }
+
+  return item;
+}
+
+static void rna_NodeGeometrySimulationOutput_items_remove(
+    ID *id, bNode *node, Main *bmain, ReportList *reports, NodeSimulationItem *item)
+{
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+  if (!NOD_geometry_simulation_output_contains_item(sim, item)) {
+    BKE_reportf(reports, RPT_ERROR, "Unable to locate item '%s' in node", item->name);
+  }
+  else {
+    NOD_geometry_simulation_output_remove_item(sim, item);
+
+    bNodeTree *ntree = (bNodeTree *)id;
+    BKE_ntree_update_tag_node_property(ntree, node);
+    ED_node_tree_propagate_change(NULL, bmain, ntree);
+    WM_main_add_notifier(NC_NODE | NA_EDITED, ntree);
+  }
+}
+
+static void rna_NodeGeometrySimulationOutput_items_clear(ID *id, bNode *node, Main *bmain)
+{
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+  NOD_geometry_simulation_output_clear_items(sim);
+
+  bNodeTree *ntree = (bNodeTree *)id;
+  BKE_ntree_update_tag_node_property(ntree, node);
+  ED_node_tree_propagate_change(NULL, bmain, ntree);
+  WM_main_add_notifier(NC_NODE | NA_EDITED, ntree);
+}
+
+static void rna_NodeGeometrySimulationOutput_items_move(
+    ID *id, bNode *node, Main *bmain, int from_index, int to_index)
+{
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+
+  if (from_index < 0 || from_index >= sim->items_num || to_index < 0 || to_index >= sim->items_num)
+  {
+    return;
+  }
+
+  NOD_geometry_simulation_output_move_item(sim, from_index, to_index);
+
+  bNodeTree *ntree = (bNodeTree *)id;
+  BKE_ntree_update_tag_node_property(ntree, node);
+  ED_node_tree_propagate_change(NULL, bmain, ntree);
+  WM_main_add_notifier(NC_NODE | NA_EDITED, ntree);
+}
+
+static PointerRNA rna_NodeGeometrySimulationOutput_active_item_get(PointerRNA *ptr)
+{
+  bNode *node = (bNode *)ptr->data;
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+  NodeSimulationItem *item = NOD_geometry_simulation_output_get_active_item(sim);
+  PointerRNA r_ptr;
+  RNA_pointer_create(ptr->owner_id, &RNA_SimulationStateItem, item, &r_ptr);
+  return r_ptr;
+}
+
+static void rna_NodeGeometrySimulationOutput_active_item_set(PointerRNA *ptr,
+                                                             PointerRNA value,
+                                                             ReportList *UNUSED(reports))
+{
+  bNode *node = (bNode *)ptr->data;
+  NodeGeometrySimulationOutput *sim = (NodeGeometrySimulationOutput *)node->storage;
+  NOD_geometry_simulation_output_set_active_item(sim, (NodeSimulationItem *)value.data);
+}
+
 /* ******** Node Socket Types ******** */
 
 static PointerRNA rna_NodeOutputFile_slot_layer_get(CollectionPropertyIterator *iter)
@@ -9739,6 +9902,150 @@ static void def_geo_set_curve_normal(StructRNA *srna)
   RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_Node_update");
 }
 
+static void def_geo_simulation_input(StructRNA *srna)
+{
+  PropertyRNA *prop;
+  FunctionRNA *func;
+  PropertyRNA *parm;
+
+  RNA_def_struct_sdna_from(srna, "NodeGeometrySimulationInput", "storage");
+
+  prop = RNA_def_property(srna, "paired_output", PROP_POINTER, PROP_NONE);
+  RNA_def_property_struct_type(prop, "Node");
+  RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+  RNA_def_property_pointer_funcs(
+      prop, "rna_NodeGeometrySimulationInput_paired_output_get", NULL, NULL, NULL);
+  RNA_def_property_ui_text(
+      prop, "Paired Output", "Simulation output node that this input node is paired with");
+
+  func = RNA_def_function(
+      srna, "pair_with_output", "rna_GeometryNodeSimulationInput_pair_with_output");
+  RNA_def_function_ui_description(func, "Pair a simulation input node with an output node.");
+  RNA_def_function_flag(func, FUNC_USE_SELF_ID | FUNC_USE_REPORTS | FUNC_USE_CONTEXT);
+  parm = RNA_def_pointer(
+      func, "output_node", "GeometryNode", "Output Node", "Simulation output node to pair with");
+  RNA_def_parameter_flags(parm, 0, PARM_REQUIRED);
+  /* return value */
+  parm = RNA_def_boolean(
+      func, "result", false, "Result", "True if pairing the node was successful");
+  RNA_def_function_return(func, parm);
+}
+
+static void rna_def_simulation_state_item(BlenderRNA *brna)
+{
+  PropertyRNA *prop;
+
+  StructRNA *srna = RNA_def_struct(brna, "SimulationStateItem", NULL);
+  RNA_def_struct_ui_text(srna, "Simulation Item", "");
+  RNA_def_struct_sdna(srna, "NodeSimulationItem");
+
+  prop = RNA_def_property(srna, "name", PROP_STRING, PROP_NONE);
+  RNA_def_property_string_funcs(prop, NULL, NULL, "rna_SimulationStateItem_name_set");
+  RNA_def_property_ui_text(prop, "Name", "");
+  RNA_def_struct_name_property(srna, prop);
+  RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_SimulationStateItem_update");
+
+  prop = RNA_def_property(srna, "socket_type", PROP_ENUM, PROP_NONE);
+  RNA_def_property_enum_items(prop, node_socket_data_type_items);
+  RNA_def_property_enum_funcs(prop, NULL, NULL, "rna_SimulationStateItem_socket_type_itemf");
+  RNA_def_property_ui_text(prop, "Socket Type", "");
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_SimulationStateItem_update");
+
+  prop = RNA_def_property(srna, "attribute_domain", PROP_ENUM, PROP_NONE);
+  RNA_def_property_enum_items(prop, rna_enum_attribute_domain_items);
+  RNA_def_property_ui_text(
+      prop,
+      "Attribute Domain",
+      "Attribute domain where the attribute domain is stored in the simulation state");
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  RNA_def_property_update(prop, NC_NODE | NA_EDITED, "rna_SimulationStateItem_update");
+
+  prop = RNA_def_property(srna, "color", PROP_FLOAT, PROP_COLOR_GAMMA);
+  RNA_def_property_array(prop, 4);
+  RNA_def_property_float_funcs(prop, "rna_SimulationStateItem_color_get", NULL, NULL);
+  RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+  RNA_def_property_ui_text(
+      prop, "Color", "Color of the corresponding socket type in the node editor");
+}
+
+static void rna_def_geo_simulation_output_items(BlenderRNA *brna)
+{
+  StructRNA *srna;
+  PropertyRNA *parm;
+  FunctionRNA *func;
+
+  srna = RNA_def_struct(brna, "NodeGeometrySimulationOutputItems", NULL);
+  RNA_def_struct_sdna(srna, "bNode");
+  RNA_def_struct_ui_text(srna, "Items", "Collection of simulation items");
+
+  func = RNA_def_function(srna, "new", "rna_NodeGeometrySimulationOutput_items_new");
+  RNA_def_function_ui_description(func, "Add a item to this simulation zone");
+  RNA_def_function_flag(func, FUNC_USE_SELF_ID | FUNC_USE_MAIN | FUNC_USE_REPORTS);
+  parm = RNA_def_enum(func,
+                      "socket_type",
+                      node_socket_data_type_items,
+                      SOCK_GEOMETRY,
+                      "Socket Type",
+                      "Socket type of the item");
+  RNA_def_parameter_flags(parm, 0, PARM_REQUIRED);
+  parm = RNA_def_string(func, "name", NULL, MAX_NAME, "Name", "");
+  RNA_def_parameter_flags(parm, 0, PARM_REQUIRED);
+  /* return value */
+  parm = RNA_def_pointer(func, "item", "SimulationStateItem", "Item", "New item");
+  RNA_def_function_return(func, parm);
+
+  func = RNA_def_function(srna, "remove", "rna_NodeGeometrySimulationOutput_items_remove");
+  RNA_def_function_ui_description(func, "Remove an item from this simulation zone");
+  RNA_def_function_flag(func, FUNC_USE_SELF_ID | FUNC_USE_MAIN | FUNC_USE_REPORTS);
+  parm = RNA_def_pointer(func, "item", "SimulationStateItem", "Item", "The item to remove");
+  RNA_def_parameter_flags(parm, PROP_NEVER_NULL, PARM_REQUIRED);
+
+  func = RNA_def_function(srna, "clear", "rna_NodeGeometrySimulationOutput_items_clear");
+  RNA_def_function_ui_description(func, "Remove all items from this simulation zone");
+  RNA_def_function_flag(func, FUNC_USE_SELF_ID | FUNC_USE_MAIN);
+
+  func = RNA_def_function(srna, "move", "rna_NodeGeometrySimulationOutput_items_move");
+  RNA_def_function_ui_description(func, "Move an item to another position");
+  RNA_def_function_flag(func, FUNC_USE_SELF_ID | FUNC_USE_MAIN);
+  parm = RNA_def_int(
+      func, "from_index", -1, 0, INT_MAX, "From Index", "Index of the item to move", 0, 10000);
+  RNA_def_parameter_flags(parm, 0, PARM_REQUIRED);
+  parm = RNA_def_int(
+      func, "to_index", -1, 0, INT_MAX, "To Index", "Target index for the item", 0, 10000);
+  RNA_def_parameter_flags(parm, 0, PARM_REQUIRED);
+}
+
+static void def_geo_simulation_output(StructRNA *srna)
+{
+  PropertyRNA *prop;
+
+  RNA_def_struct_sdna_from(srna, "NodeGeometrySimulationOutput", "storage");
+
+  prop = RNA_def_property(srna, "state_items", PROP_COLLECTION, PROP_NONE);
+  RNA_def_property_collection_sdna(prop, NULL, "items", "items_num");
+  RNA_def_property_struct_type(prop, "SimulationStateItem");
+  RNA_def_property_ui_text(prop, "Items", "");
+  RNA_def_property_srna(prop, "NodeGeometrySimulationOutputItems");
+
+  prop = RNA_def_property(srna, "active_index", PROP_INT, PROP_UNSIGNED);
+  RNA_def_property_int_sdna(prop, NULL, "active_index");
+  RNA_def_property_ui_text(prop, "Active Item Index", "Index of the active item");
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  RNA_def_property_update(prop, NC_NODE, NULL);
+
+  prop = RNA_def_property(srna, "active_item", PROP_POINTER, PROP_NONE);
+  RNA_def_property_struct_type(prop, "SimulationStateItem");
+  RNA_def_property_pointer_funcs(prop,
+                                 "rna_NodeGeometrySimulationOutput_active_item_get",
+                                 "rna_NodeGeometrySimulationOutput_active_item_set",
+                                 NULL,
+                                 NULL);
+  RNA_def_property_flag(prop, PROP_EDITABLE);
+  RNA_def_property_ui_text(prop, "Active Item Index", "Index of the active item");
+  RNA_def_property_update(prop, NC_NODE, NULL);
+}
+
 static void def_geo_curve_handle_type_selection(StructRNA *srna)
 {
   PropertyRNA *prop;
@@ -13098,6 +13405,7 @@ void RNA_def_nodetree(BlenderRNA *brna)
   rna_def_compositor_node(brna);
   rna_def_texture_node(brna);
   rna_def_geometry_node(brna);
+  rna_def_simulation_state_item(brna);
   rna_def_function_node(brna);
 
   rna_def_nodetree(brna);
@@ -13159,6 +13467,7 @@ void RNA_def_nodetree(BlenderRNA *brna)
   /* special socket types */
   rna_def_cmp_output_file_slot_file(brna);
   rna_def_cmp_output_file_slot_layer(brna);
+  rna_def_geo_simulation_output_items(brna);
 
   rna_def_node_instance_hash(brna);
 }

source/blender/makesrna/intern/rna_userdef.c

@@ -3057,6 +3057,12 @@ static void rna_def_userdef_theme_space_node(BlenderRNA *brna)
   RNA_def_property_array(prop, 3);
   RNA_def_property_ui_text(prop, "Attribute Node", "");
   RNA_def_property_update(prop, 0, "rna_userdef_theme_update");
+
+  prop = RNA_def_property(srna, "simulation_zone", PROP_FLOAT, PROP_COLOR_GAMMA);
+  RNA_def_property_float_sdna(prop, NULL, "node_zone_simulation");
+  RNA_def_property_array(prop, 4);
+  RNA_def_property_ui_text(prop, "Simulation Zone", "");
+  RNA_def_property_update(prop, 0, "rna_userdef_theme_update");
 }
 
 static void rna_def_userdef_theme_space_buts(BlenderRNA *brna)

source/blender/modifiers/intern/MOD_nodes.cc

@@ -54,6 +54,8 @@
 #include "BKE_pointcloud.h"
 #include "BKE_screen.h"
 #include "BKE_simulation.h"
+#include "BKE_simulation_state.hh"
+#include "BKE_simulation_state_serialize.hh"
 #include "BKE_workspace.h"
 
 #include "BLO_read_write.h"
@@ -304,6 +306,9 @@ static bool check_tree_for_time_node(const bNodeTree &tree, Set<const bNodeTree
   if (!tree.nodes_by_type("GeometryNodeInputSceneTime").is_empty()) {
     return true;
   }
+  if (!tree.nodes_by_type("GeometryNodeSimulationInput").is_empty()) {
+    return true;
+  }
   for (const bNode *node : tree.group_nodes()) {
     if (const bNodeTree *sub_tree = reinterpret_cast<const bNodeTree *>(node->id)) {
       if (check_tree_for_time_node(*sub_tree, checked_groups)) {
@@ -1125,6 +1130,96 @@ static void store_output_attributes(GeometrySet &geometry,
   store_computed_output_attributes(geometry, attributes_to_store);
 }
 
+static void prepare_simulation_states_for_evaluation(const NodesModifierData &nmd,
+                                                     NodesModifierData &nmd_orig,
+                                                     const ModifierEvalContext &ctx,
+                                                     nodes::GeoNodesModifierData &exec_data)
+{
+  const Main *bmain = DEG_get_bmain(ctx.depsgraph);
+  const SubFrame current_frame = DEG_get_ctime(ctx.depsgraph);
+  const Scene *scene = DEG_get_input_scene(ctx.depsgraph);
+  const SubFrame start_frame = scene->r.sfra;
+  const bool is_start_frame = current_frame == start_frame;
+
+  if (DEG_is_active(ctx.depsgraph)) {
+    if (nmd_orig.simulation_cache == nullptr) {
+      nmd_orig.simulation_cache = MEM_new<bke::sim::ModifierSimulationCache>(__func__);
+    }
+
+    {
+      /* Try to use baked data. */
+      const StringRefNull bmain_path = BKE_main_blendfile_path(bmain);
+      if (nmd_orig.simulation_cache->cache_state() != bke::sim::CacheState::Baked &&
+          !bmain_path.is_empty())
+      {
+        nmd_orig.simulation_cache->try_discover_bake(
+            bke::sim::get_meta_directory(*bmain, *ctx.object, nmd.modifier),
+            bke::sim::get_bdata_directory(*bmain, *ctx.object, nmd.modifier));
+      }
+    }
+
+    {
+      /* Reset cached data if necessary. */
+      const bke::sim::StatesAroundFrame sim_states =
+          nmd_orig.simulation_cache->get_states_around_frame(current_frame);
+      if (nmd_orig.simulation_cache->cache_state() == bke::sim::CacheState::Invalid &&
+          (current_frame == start_frame ||
+           (sim_states.current == nullptr && sim_states.prev == nullptr &&
+            sim_states.next != nullptr)))
+      {
+        nmd_orig.simulation_cache->reset();
+      }
+    }
+    /* Decide if a new simulation state should be created in this evaluation. */
+    const bke::sim::StatesAroundFrame sim_states =
+        nmd_orig.simulation_cache->get_states_around_frame(current_frame);
+    if (nmd_orig.simulation_cache->cache_state() != bke::sim::CacheState::Baked) {
+      if (sim_states.current == nullptr) {
+        if (is_start_frame || !nmd_orig.simulation_cache->has_states()) {
+          bke::sim::ModifierSimulationState &current_sim_state =
+              nmd_orig.simulation_cache->get_state_at_frame_for_write(current_frame);
+          exec_data.current_simulation_state_for_write = &current_sim_state;
+          exec_data.simulation_time_delta = 0.0f;
+          if (!is_start_frame) {
+            /* When starting a new simulation at another frame than the start frame, it can't match
+             * what would be baked, so invalidate it immediately. */
+            nmd_orig.simulation_cache->invalidate();
+          }
+        }
+        else if (sim_states.prev != nullptr && sim_states.next == nullptr) {
+          const float delta_frames = float(current_frame) - float(sim_states.prev->frame);
+          if (delta_frames <= 1.0f) {
+            bke::sim::ModifierSimulationState &current_sim_state =
+                nmd_orig.simulation_cache->get_state_at_frame_for_write(current_frame);
+            exec_data.current_simulation_state_for_write = &current_sim_state;
+            const float delta_seconds = delta_frames / FPS;
+            exec_data.simulation_time_delta = delta_seconds;
+          }
+        }
+      }
+    }
+  }
+
+  /* Load read-only states to give nodes access to cached data. */
+  const bke::sim::StatesAroundFrame sim_states =
+      nmd_orig.simulation_cache->get_states_around_frame(current_frame);
+  if (sim_states.current) {
+    sim_states.current->state.ensure_bake_loaded();
+    exec_data.current_simulation_state = &sim_states.current->state;
+  }
+  if (sim_states.prev) {
+    sim_states.prev->state.ensure_bake_loaded();
+    exec_data.prev_simulation_state = &sim_states.prev->state;
+    if (sim_states.next) {
+      sim_states.next->state.ensure_bake_loaded();
+      exec_data.next_simulation_state = &sim_states.next->state;
+      exec_data.simulation_state_mix_factor =
+          (float(current_frame) - float(sim_states.prev->frame)) /
+          (float(sim_states.next->frame) - float(sim_states.prev->frame));
+    }
+  }
+}
+
 /**
  * Evaluate a node group to compute the output geometry.
  */
@@ -1135,6 +1230,9 @@ static GeometrySet compute_geometry(const bNodeTree &btree,
                                     NodesModifierData *nmd,
                                     const ModifierEvalContext *ctx)
 {
+  NodesModifierData *nmd_orig = reinterpret_cast<NodesModifierData *>(
+      BKE_modifier_get_original(ctx->object, &nmd->modifier));
+
   const nodes::GeometryNodeLazyFunctionGraphMapping &mapping = lf_graph_info.mapping;
 
   Vector<const lf::OutputSocket *> graph_inputs = mapping.group_input_sockets;
@@ -1160,6 +1258,8 @@ static GeometrySet compute_geometry(const bNodeTree &btree,
   geo_nodes_modifier_data.self_object = ctx->object;
   auto eval_log = std::make_unique<geo_log::GeoModifierLog>();
 
+  prepare_simulation_states_for_evaluation(*nmd, *nmd_orig, *ctx, geo_nodes_modifier_data);
+
   Set<ComputeContextHash> socket_log_contexts;
   if (logging_enabled(ctx)) {
     geo_nodes_modifier_data.eval_log = eval_log.get();
@@ -1240,8 +1340,6 @@ static GeometrySet compute_geometry(const bNodeTree &btree,
   }
 
   if (logging_enabled(ctx)) {
-    NodesModifierData *nmd_orig = reinterpret_cast<NodesModifierData *>(
-        BKE_modifier_get_original(ctx->object, &nmd->modifier));
     delete static_cast<geo_log::GeoModifierLog *>(nmd_orig->runtime_eval_log);
     nmd_orig->runtime_eval_log = eval_log.release();
   }
@@ -1928,6 +2026,7 @@ static void blendRead(BlendDataReader *reader, ModifierData *md)
     IDP_BlendDataRead(reader, &nmd->settings.properties);
   }
   nmd->runtime_eval_log = nullptr;
+  nmd->simulation_cache = nullptr;
 }
 
 static void copyData(const ModifierData *md, ModifierData *target, const int flag)
@@ -1938,6 +2037,7 @@ static void copyData(const ModifierData *md, ModifierData *target, const int fla
   BKE_modifier_copydata_generic(md, target, flag);
 
   tnmd->runtime_eval_log = nullptr;
+  tnmd->simulation_cache = nullptr;
 
   if (nmd->settings.properties != nullptr) {
     tnmd->settings.properties = IDP_CopyProperty_ex(nmd->settings.properties, flag);
@@ -1952,6 +2052,8 @@ static void freeData(ModifierData *md)
     nmd->settings.properties = nullptr;
   }
 
+  MEM_delete(nmd->simulation_cache);
+
   clear_runtime_data(nmd);
 }
 

source/blender/nodes/NOD_geometry.h

@@ -12,6 +12,74 @@ extern struct bNodeTreeType *ntreeType_Geometry;
 
 void register_node_type_geo_custom_group(bNodeType *ntype);
 
+/* -------------------------------------------------------------------- */
+/** \name Simulation Input Node
+ * \{ */
+
+struct bNode *NOD_geometry_simulation_input_get_paired_output(
+    struct bNodeTree *node_tree, const struct bNode *simulation_input_node);
+
+/**
+ * Pair a simulation input node with an output node.
+ * \return True if pairing the node was successful.
+ */
+bool NOD_geometry_simulation_input_pair_with_output(const struct bNodeTree *node_tree,
+                                                    struct bNode *simulation_input_node,
+                                                    const struct bNode *simulation_output_node);
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Simulation Output Node
+ * \{ */
+
+bool NOD_geometry_simulation_output_item_socket_type_supported(eNodeSocketDatatype socket_type);
+
+/**
+ * Set a unique item name.
+ * \return True if the unique name differs from the original name.
+ */
+bool NOD_geometry_simulation_output_item_set_unique_name(struct NodeGeometrySimulationOutput *sim,
+                                                         struct NodeSimulationItem *item,
+                                                         const char *name,
+                                                         const char *defname);
+
+/**
+ * Find the node owning this simulation state item.
+ */
+bNode *NOD_geometry_simulation_output_find_node_by_item(struct bNodeTree *ntree,
+                                                        const struct NodeSimulationItem *item);
+
+bool NOD_geometry_simulation_output_contains_item(struct NodeGeometrySimulationOutput *sim,
+                                                  const struct NodeSimulationItem *item);
+struct NodeSimulationItem *NOD_geometry_simulation_output_get_active_item(
+    struct NodeGeometrySimulationOutput *sim);
+void NOD_geometry_simulation_output_set_active_item(struct NodeGeometrySimulationOutput *sim,
+                                                    struct NodeSimulationItem *item);
+struct NodeSimulationItem *NOD_geometry_simulation_output_find_item(
+    struct NodeGeometrySimulationOutput *sim, const char *name);
+struct NodeSimulationItem *NOD_geometry_simulation_output_add_item(
+    struct NodeGeometrySimulationOutput *sim, short socket_type, const char *name);
+struct NodeSimulationItem *NOD_geometry_simulation_output_insert_item(
+    struct NodeGeometrySimulationOutput *sim, short socket_type, const char *name, int index);
+struct NodeSimulationItem *NOD_geometry_simulation_output_add_item_from_socket(
+    struct NodeGeometrySimulationOutput *sim,
+    const struct bNode *from_node,
+    const struct bNodeSocket *from_sock);
+struct NodeSimulationItem *NOD_geometry_simulation_output_insert_item_from_socket(
+    struct NodeGeometrySimulationOutput *sim,
+    const struct bNode *from_node,
+    const struct bNodeSocket *from_sock,
+    int index);
+void NOD_geometry_simulation_output_remove_item(struct NodeGeometrySimulationOutput *sim,
+                                                struct NodeSimulationItem *item);
+void NOD_geometry_simulation_output_clear_items(struct NodeGeometrySimulationOutput *sim);
+void NOD_geometry_simulation_output_move_item(struct NodeGeometrySimulationOutput *sim,
+                                              int from_index,
+                                              int to_index);
+
+/** \} */
+
 #ifdef __cplusplus
 }
 #endif

source/blender/nodes/NOD_geometry_nodes_lazy_function.hh

@@ -26,6 +26,8 @@
 
 #include "BLI_compute_context.hh"
 
+#include "BKE_simulation_state.hh"
+
 struct Object;
 struct Depsgraph;
 
@@ -44,6 +46,16 @@ struct GeoNodesModifierData {
   Depsgraph *depsgraph = nullptr;
   /** Optional logger. */
   geo_eval_log::GeoModifierLog *eval_log = nullptr;
+
+  /** Read-only simulation states around the current frame. */
+  const bke::sim::ModifierSimulationState *current_simulation_state = nullptr;
+  const bke::sim::ModifierSimulationState *prev_simulation_state = nullptr;
+  const bke::sim::ModifierSimulationState *next_simulation_state = nullptr;
+  float simulation_state_mix_factor = 0.0f;
+  /** Used when the evaluation should create a new simulation state. */
+  bke::sim::ModifierSimulationState *current_simulation_state_for_write = nullptr;
+  float simulation_time_delta = 0.0f;
+
   /**
    * Some nodes should be executed even when their output is not used (e.g. active viewer nodes and
    * the node groups they are contained in).
@@ -149,6 +161,7 @@ struct GeometryNodeLazyFunctionGraphMapping {
    */
   Map<const bNode *, const lf::FunctionNode *> group_node_map;
   Map<const bNode *, const lf::FunctionNode *> viewer_node_map;
+  Map<const bNode *, const lf::FunctionNode *> sim_output_node_map;
 
   /* Indexed by #bNodeSocket::index_in_all_outputs. */
   Array<int> lf_input_index_for_output_bsocket_usage;
@@ -228,8 +241,34 @@ class GeometryNodesLazyFunctionLogger : public fn::lazy_function::GraphExecutor:
                                const lf::Context &context) const override;
 };
 
+std::unique_ptr<LazyFunction> get_simulation_output_lazy_function(
+    const bNode &node, GeometryNodesLazyFunctionGraphInfo &own_lf_graph_info);
+std::unique_ptr<LazyFunction> get_simulation_input_lazy_function(
+    const bNodeTree &node_tree,
+    const bNode &node,
+    GeometryNodesLazyFunctionGraphInfo &own_lf_graph_info);
 std::unique_ptr<LazyFunction> get_switch_node_lazy_function(const bNode &node);
 
+bke::sim::SimulationZoneID get_simulation_zone_id(const ComputeContext &context,
+                                                  const int output_node_id);
+
+/**
+ * An anonymous attribute created by a node.
+ */
+class NodeAnonymousAttributeID : public bke::AnonymousAttributeID {
+  std::string long_name_;
+  std::string socket_name_;
+
+ public:
+  NodeAnonymousAttributeID(const Object &object,
+                           const ComputeContext &compute_context,
+                           const bNode &bnode,
+                           const StringRef identifier,
+                           const StringRef name);
+
+  std::string user_name() const override;
+};
+
 /**
  * Tells the lazy-function graph evaluator which nodes have side effects based on the current
  * context. For example, the same viewer node can have side effects in one context, but not in

source/blender/nodes/NOD_static_types.h

@@ -419,6 +419,8 @@ DefNode(GeometryNode, GEO_NODE_SET_POSITION, 0, "SET_POSITION", SetPosition, "Se
 DefNode(GeometryNode, GEO_NODE_SET_SHADE_SMOOTH, 0, "SET_SHADE_SMOOTH", SetShadeSmooth, "Set Shade Smooth", "Control the smoothness of mesh normals around each face by changing the \"shade smooth\" attribute")
 DefNode(GeometryNode, GEO_NODE_SET_SPLINE_CYCLIC, 0, "SET_SPLINE_CYCLIC", SetSplineCyclic, "Set Spline Cyclic", "Control whether each spline loops back on itself by changing the \"cyclic\" attribute")
 DefNode(GeometryNode, GEO_NODE_SET_SPLINE_RESOLUTION, 0, "SET_SPLINE_RESOLUTION", SetSplineResolution, "Set Spline Resolution", "Control how many evaluated points should be generated on every curve segment")
+DefNode(GeometryNode, GEO_NODE_SIMULATION_INPUT, def_geo_simulation_input, "SIMULATION_INPUT", SimulationInput, "Simulation Input", "")
+DefNode(GeometryNode, GEO_NODE_SIMULATION_OUTPUT, def_geo_simulation_output, "SIMULATION_OUTPUT", SimulationOutput, "Simulation Output", "")
 DefNode(GeometryNode, GEO_NODE_SPLIT_EDGES, 0, "SPLIT_EDGES", SplitEdges, "Split Edges", "Duplicate mesh edges and break connections with the surrounding faces")
 DefNode(GeometryNode, GEO_NODE_STORE_NAMED_ATTRIBUTE, def_geo_store_named_attribute, "STORE_NAMED_ATTRIBUTE", StoreNamedAttribute, "Store Named Attribute", "Store the result of a field on a geometry as an attribute with the specified name")
 DefNode(GeometryNode, GEO_NODE_STRING_JOIN, 0, "STRING_JOIN", StringJoin, "Join Strings", "Combine any number of input strings")

source/blender/nodes/geometry/CMakeLists.txt

@@ -170,6 +170,8 @@ set(SRC
   nodes/node_geo_set_shade_smooth.cc
   nodes/node_geo_set_spline_cyclic.cc
   nodes/node_geo_set_spline_resolution.cc
+  nodes/node_geo_simulation_input.cc
+  nodes/node_geo_simulation_output.cc
   nodes/node_geo_store_named_attribute.cc
   nodes/node_geo_string_join.cc
   nodes/node_geo_string_to_curves.cc

source/blender/nodes/geometry/node_geometry_register.cc

@@ -154,6 +154,8 @@ void register_geometry_nodes()
   register_node_type_geo_set_shade_smooth();
   register_node_type_geo_set_spline_cyclic();
   register_node_type_geo_set_spline_resolution();
+  register_node_type_geo_simulation_input();
+  register_node_type_geo_simulation_output();
   register_node_type_geo_store_named_attribute();
   register_node_type_geo_string_join();
   register_node_type_geo_string_to_curves();

source/blender/nodes/geometry/node_geometry_register.hh

@@ -152,6 +152,8 @@ void register_node_type_geo_set_position();
 void register_node_type_geo_set_shade_smooth();
 void register_node_type_geo_set_spline_cyclic();
 void register_node_type_geo_set_spline_resolution();
+void register_node_type_geo_simulation_input();
+void register_node_type_geo_simulation_output();
 void register_node_type_geo_store_named_attribute();
 void register_node_type_geo_string_join();
 void register_node_type_geo_string_to_curves();

source/blender/nodes/geometry/node_geometry_util.hh

@@ -133,4 +133,20 @@ class FieldAtIndexInput final : public bke::GeometryFieldInput {
   }
 };
 
+std::string socket_identifier_for_simulation_item(const NodeSimulationItem &item);
+
+void socket_declarations_for_simulation_items(Span<NodeSimulationItem> items,
+                                              NodeDeclaration &r_declaration);
+const CPPType &get_simulation_item_cpp_type(eNodeSocketDatatype socket_type);
+const CPPType &get_simulation_item_cpp_type(const NodeSimulationItem &item);
+void values_to_simulation_state(const Span<NodeSimulationItem> node_simulation_items,
+                                const Span<void *> input_values,
+                                bke::sim::SimulationZoneState &r_zone_state);
+void simulation_state_to_values(const Span<NodeSimulationItem> node_simulation_items,
+                                const bke::sim::SimulationZoneState &zone_state,
+                                const Object &self_object,
+                                const ComputeContext &compute_context,
+                                const bNode &sim_output_node,
+                                Span<void *> r_output_values);
+
 }  // namespace blender::nodes

source/blender/nodes/geometry/nodes/node_geo_simulation_input.cc

@@ -0,0 +1,253 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BKE_compute_contexts.hh"
+#include "BKE_scene.h"
+
+#include "DEG_depsgraph_query.h"
+
+#include "UI_interface.h"
+#include "UI_resources.h"
+
+#include "NOD_geometry.h"
+#include "NOD_socket.h"
+
+#include "node_geometry_util.hh"
+
+namespace blender::nodes::node_geo_simulation_input_cc {
+
+NODE_STORAGE_FUNCS(NodeGeometrySimulationInput);
+
+class LazyFunctionForSimulationInputNode final : public LazyFunction {
+  const bNode &node_;
+  int32_t output_node_id_;
+  Span<NodeSimulationItem> simulation_items_;
+
+ public:
+  LazyFunctionForSimulationInputNode(const bNodeTree &node_tree,
+                                     const bNode &node,
+                                     GeometryNodesLazyFunctionGraphInfo &own_lf_graph_info)
+      : node_(node)
+  {
+    debug_name_ = "Simulation Input";
+    output_node_id_ = node_storage(node).output_node_id;
+    const bNode &output_node = *node_tree.node_by_id(output_node_id_);
+    const NodeGeometrySimulationOutput &storage = *static_cast<NodeGeometrySimulationOutput *>(
+        output_node.storage);
+    simulation_items_ = {storage.items, storage.items_num};
+
+    MutableSpan<int> lf_index_by_bsocket = own_lf_graph_info.mapping.lf_index_by_bsocket;
+    lf_index_by_bsocket[node.output_socket(0).index_in_tree()] = outputs_.append_and_get_index_as(
+        "Delta Time", CPPType::get<ValueOrField<float>>());
+
+    for (const int i : simulation_items_.index_range()) {
+      const NodeSimulationItem &item = simulation_items_[i];
+      const bNodeSocket &input_bsocket = node.input_socket(i);
+      const bNodeSocket &output_bsocket = node.output_socket(i + 1);
+
+      const CPPType &type = get_simulation_item_cpp_type(item);
+
+      lf_index_by_bsocket[input_bsocket.index_in_tree()] = inputs_.append_and_get_index_as(
+          item.name, type, lf::ValueUsage::Maybe);
+      lf_index_by_bsocket[output_bsocket.index_in_tree()] = outputs_.append_and_get_index_as(
+          item.name, type);
+    }
+  }
+
+  void execute_impl(lf::Params &params, const lf::Context &context) const final
+  {
+    const GeoNodesLFUserData &user_data = *static_cast<const GeoNodesLFUserData *>(
+        context.user_data);
+    const GeoNodesModifierData &modifier_data = *user_data.modifier_data;
+    if (modifier_data.current_simulation_state == nullptr) {
+      params.set_default_remaining_outputs();
+      return;
+    }
+
+    if (!params.output_was_set(0)) {
+      const float delta_time = modifier_data.simulation_time_delta;
+      params.set_output(0, fn::ValueOrField<float>(delta_time));
+    }
+
+    const bke::sim::SimulationZoneID zone_id = get_simulation_zone_id(*user_data.compute_context,
+                                                                      output_node_id_);
+
+    const bke::sim::SimulationZoneState *prev_zone_state =
+        modifier_data.prev_simulation_state == nullptr ?
+            nullptr :
+            modifier_data.prev_simulation_state->get_zone_state(zone_id);
+
+    std::optional<bke::sim::SimulationZoneState> initial_prev_zone_state;
+    if (prev_zone_state == nullptr) {
+      Array<void *> input_values(simulation_items_.size(), nullptr);
+      for (const int i : simulation_items_.index_range()) {
+        input_values[i] = params.try_get_input_data_ptr_or_request(i);
+      }
+      if (input_values.as_span().contains(nullptr)) {
+        /* Wait until all inputs are available. */
+        return;
+      }
+
+      initial_prev_zone_state.emplace();
+      values_to_simulation_state(simulation_items_, input_values, *initial_prev_zone_state);
+      prev_zone_state = &*initial_prev_zone_state;
+    }
+
+    Array<void *> output_values(simulation_items_.size());
+    for (const int i : simulation_items_.index_range()) {
+      output_values[i] = params.get_output_data_ptr(i + 1);
+    }
+    simulation_state_to_values(simulation_items_,
+                               *prev_zone_state,
+                               *modifier_data.self_object,
+                               *user_data.compute_context,
+                               node_,
+                               output_values);
+    for (const int i : simulation_items_.index_range()) {
+      params.output_set(i + 1);
+    }
+  }
+};
+
+}  // namespace blender::nodes::node_geo_simulation_input_cc
+
+namespace blender::nodes {
+
+std::unique_ptr<LazyFunction> get_simulation_input_lazy_function(
+    const bNodeTree &node_tree,
+    const bNode &node,
+    GeometryNodesLazyFunctionGraphInfo &own_lf_graph_info)
+{
+  namespace file_ns = blender::nodes::node_geo_simulation_input_cc;
+  BLI_assert(node.type == GEO_NODE_SIMULATION_INPUT);
+  return std::make_unique<file_ns::LazyFunctionForSimulationInputNode>(
+      node_tree, node, own_lf_graph_info);
+}
+
+}  // namespace blender::nodes
+
+namespace blender::nodes::node_geo_simulation_input_cc {
+
+static void node_declare_dynamic(const bNodeTree &node_tree,
+                                 const bNode &node,
+                                 NodeDeclaration &r_declaration)
+{
+  const bNode *output_node = node_tree.node_by_id(node_storage(node).output_node_id);
+  if (!output_node) {
+    return;
+  }
+
+  std::unique_ptr<decl::Float> delta_time = std::make_unique<decl::Float>();
+  delta_time->identifier = "Delta Time";
+  delta_time->name = DATA_("Delta Time");
+  delta_time->in_out = SOCK_OUT;
+  r_declaration.outputs.append(std::move(delta_time));
+
+  const NodeGeometrySimulationOutput &storage = *static_cast<const NodeGeometrySimulationOutput *>(
+      output_node->storage);
+  socket_declarations_for_simulation_items({storage.items, storage.items_num}, r_declaration);
+}
+
+static void node_init(bNodeTree * /*tree*/, bNode *node)
+{
+  NodeGeometrySimulationInput *data = MEM_cnew<NodeGeometrySimulationInput>(__func__);
+  /* Needs to be initialized for the node to work. */
+  data->output_node_id = 0;
+  node->storage = data;
+}
+
+static bool node_insert_link(bNodeTree *ntree, bNode *node, bNodeLink *link)
+{
+  const bNode *output_node = ntree->node_by_id(node_storage(*node).output_node_id);
+  if (!output_node) {
+    return true;
+  }
+
+  NodeGeometrySimulationOutput &storage = *static_cast<NodeGeometrySimulationOutput *>(
+      output_node->storage);
+
+  if (link->tonode == node) {
+    if (link->tosock->identifier == StringRef("__extend__")) {
+      if (const NodeSimulationItem *item = NOD_geometry_simulation_output_add_item_from_socket(
+              &storage, link->fromnode, link->fromsock))
+      {
+        update_node_declaration_and_sockets(*ntree, *node);
+        link->tosock = nodeFindSocket(
+            node, SOCK_IN, socket_identifier_for_simulation_item(*item).c_str());
+      }
+      else {
+        return false;
+      }
+    }
+  }
+  else {
+    BLI_assert(link->fromnode == node);
+    if (link->fromsock->identifier == StringRef("__extend__")) {
+      if (const NodeSimulationItem *item = NOD_geometry_simulation_output_add_item_from_socket(
+              &storage, link->tonode, link->tosock))
+      {
+        update_node_declaration_and_sockets(*ntree, *node);
+        link->fromsock = nodeFindSocket(
+            node, SOCK_OUT, socket_identifier_for_simulation_item(*item).c_str());
+      }
+      else {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+}  // namespace blender::nodes::node_geo_simulation_input_cc
+
+void register_node_type_geo_simulation_input()
+{
+  namespace file_ns = blender::nodes::node_geo_simulation_input_cc;
+
+  static bNodeType ntype;
+  geo_node_type_base(&ntype, GEO_NODE_SIMULATION_INPUT, "Simulation Input", NODE_CLASS_INTERFACE);
+  ntype.initfunc = file_ns::node_init;
+  ntype.declare_dynamic = file_ns::node_declare_dynamic;
+  ntype.insert_link = file_ns::node_insert_link;
+  ntype.gather_add_node_search_ops = nullptr;
+  ntype.gather_link_search_ops = nullptr;
+  node_type_storage(&ntype,
+                    "NodeGeometrySimulationInput",
+                    node_free_standard_storage,
+                    node_copy_standard_storage);
+  nodeRegisterType(&ntype);
+}
+
+bNode *NOD_geometry_simulation_input_get_paired_output(bNodeTree *node_tree,
+                                                       const bNode *simulation_input_node)
+{
+  namespace file_ns = blender::nodes::node_geo_simulation_input_cc;
+
+  const NodeGeometrySimulationInput &data = file_ns::node_storage(*simulation_input_node);
+  return node_tree->node_by_id(data.output_node_id);
+}
+
+bool NOD_geometry_simulation_input_pair_with_output(const bNodeTree *node_tree,
+                                                    bNode *sim_input_node,
+                                                    const bNode *sim_output_node)
+{
+  namespace file_ns = blender::nodes::node_geo_simulation_input_cc;
+
+  BLI_assert(sim_input_node->type == GEO_NODE_SIMULATION_INPUT);
+  if (sim_output_node->type != GEO_NODE_SIMULATION_OUTPUT) {
+    return false;
+  }
+
+  /* Allow only one input paired to an output. */
+  for (const bNode *other_input_node : node_tree->nodes_by_type("GeometryNodeSimulationInput")) {
+    if (other_input_node != sim_input_node) {
+      const NodeGeometrySimulationInput &other_storage = file_ns::node_storage(*other_input_node);
+      if (other_storage.output_node_id == sim_output_node->identifier) {
+        return false;
+      }
+    }
+  }
+
+  NodeGeometrySimulationInput &storage = file_ns::node_storage(*sim_input_node);
+  storage.output_node_id = sim_output_node->identifier;
+  return true;
+}

source/blender/nodes/geometry/nodes/node_geo_simulation_output.cc

@@ -0,0 +1,893 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BLI_string_utils.h"
+
+#include "BKE_compute_contexts.hh"
+#include "BKE_curves.hh"
+#include "BKE_instances.hh"
+#include "BKE_scene.h"
+
+#include "DEG_depsgraph_query.h"
+
+#include "UI_interface.h"
+
+#include "NOD_common.h"
+#include "NOD_socket.h"
+
+#include "FN_field_cpp_type.hh"
+
+#include "DNA_curves_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_pointcloud_types.h"
+
+#include "NOD_add_node_search.hh"
+
+#include "node_geometry_util.hh"
+
+namespace blender::nodes {
+
+std::string socket_identifier_for_simulation_item(const NodeSimulationItem &item)
+{
+  return "Item_" + std::to_string(item.identifier);
+}
+
+static std::unique_ptr<SocketDeclaration> socket_declaration_for_simulation_item(
+    const NodeSimulationItem &item, const eNodeSocketInOut in_out, const int index)
+{
+  const eNodeSocketDatatype socket_type = eNodeSocketDatatype(item.socket_type);
+  BLI_assert(NOD_geometry_simulation_output_item_socket_type_supported(socket_type));
+
+  std::unique_ptr<SocketDeclaration> decl;
+  switch (socket_type) {
+    case SOCK_FLOAT:
+      decl = std::make_unique<decl::Float>();
+      decl->input_field_type = InputSocketFieldType::IsSupported;
+      decl->output_field_dependency = OutputFieldDependency::ForPartiallyDependentField({index});
+      break;
+    case SOCK_VECTOR:
+      decl = std::make_unique<decl::Vector>();
+      decl->input_field_type = InputSocketFieldType::IsSupported;
+      decl->output_field_dependency = OutputFieldDependency::ForPartiallyDependentField({index});
+      break;
+    case SOCK_RGBA:
+      decl = std::make_unique<decl::Color>();
+      decl->input_field_type = InputSocketFieldType::IsSupported;
+      decl->output_field_dependency = OutputFieldDependency::ForPartiallyDependentField({index});
+      break;
+    case SOCK_BOOLEAN:
+      decl = std::make_unique<decl::Bool>();
+      decl->input_field_type = InputSocketFieldType::IsSupported;
+      decl->output_field_dependency = OutputFieldDependency::ForPartiallyDependentField({index});
+      break;
+    case SOCK_INT:
+      decl = std::make_unique<decl::Int>();
+      decl->input_field_type = InputSocketFieldType::IsSupported;
+      decl->output_field_dependency = OutputFieldDependency::ForPartiallyDependentField({index});
+      break;
+    case SOCK_STRING:
+      decl = std::make_unique<decl::String>();
+      break;
+    case SOCK_GEOMETRY:
+      decl = std::make_unique<decl::Geometry>();
+      break;
+    default:
+      BLI_assert_unreachable();
+  }
+
+  decl->name = item.name ? item.name : "";
+  decl->identifier = socket_identifier_for_simulation_item(item);
+  decl->in_out = in_out;
+  return decl;
+}
+
+void socket_declarations_for_simulation_items(const Span<NodeSimulationItem> items,
+                                              NodeDeclaration &r_declaration)
+{
+  for (const int i : items.index_range()) {
+    const NodeSimulationItem &item = items[i];
+    r_declaration.inputs.append(socket_declaration_for_simulation_item(item, SOCK_IN, i));
+    r_declaration.outputs.append(socket_declaration_for_simulation_item(item, SOCK_OUT, i));
+  }
+  r_declaration.inputs.append(decl::create_extend_declaration(SOCK_IN));
+  r_declaration.outputs.append(decl::create_extend_declaration(SOCK_OUT));
+}
+
+struct SimulationItemsUniqueNameArgs {
+  NodeGeometrySimulationOutput *sim;
+  const NodeSimulationItem *item;
+};
+
+static bool simulation_items_unique_name_check(void *arg, const char *name)
+{
+  const SimulationItemsUniqueNameArgs &args = *static_cast<const SimulationItemsUniqueNameArgs *>(
+      arg);
+  for (const NodeSimulationItem &item : args.sim->items_span()) {
+    if (&item != args.item) {
+      if (STREQ(item.name, name)) {
+        return true;
+      }
+    }
+  }
+  if (STREQ(name, "Delta Time")) {
+    return true;
+  }
+  return false;
+}
+
+const CPPType &get_simulation_item_cpp_type(const eNodeSocketDatatype socket_type)
+{
+  const char *socket_idname = nodeStaticSocketType(socket_type, 0);
+  const bNodeSocketType *typeinfo = nodeSocketTypeFind(socket_idname);
+  BLI_assert(typeinfo);
+  BLI_assert(typeinfo->geometry_nodes_cpp_type);
+  return *typeinfo->geometry_nodes_cpp_type;
+}
+
+const CPPType &get_simulation_item_cpp_type(const NodeSimulationItem &item)
+{
+  return get_simulation_item_cpp_type(eNodeSocketDatatype(item.socket_type));
+}
+
+static void cleanup_geometry_for_simulation_state(GeometrySet &main_geometry)
+{
+  main_geometry.modify_geometry_sets([&](GeometrySet &geometry) {
+    if (Mesh *mesh = geometry.get_mesh_for_write()) {
+      mesh->attributes_for_write().remove_anonymous();
+    }
+    if (Curves *curves = geometry.get_curves_for_write()) {
+      curves->geometry.wrap().attributes_for_write().remove_anonymous();
+    }
+    if (PointCloud *pointcloud = geometry.get_pointcloud_for_write()) {
+      pointcloud->attributes_for_write().remove_anonymous();
+    }
+    if (bke::Instances *instances = geometry.get_instances_for_write()) {
+      instances->attributes_for_write().remove_anonymous();
+    }
+    geometry.keep_only_during_modify({GEO_COMPONENT_TYPE_MESH,
+                                      GEO_COMPONENT_TYPE_CURVE,
+                                      GEO_COMPONENT_TYPE_POINT_CLOUD,
+                                      GEO_COMPONENT_TYPE_INSTANCES});
+  });
+}
+
+void simulation_state_to_values(const Span<NodeSimulationItem> node_simulation_items,
+                                const bke::sim::SimulationZoneState &zone_state,
+                                const Object &self_object,
+                                const ComputeContext &compute_context,
+                                const bNode &node,
+                                Span<void *> r_output_values)
+{
+  /* Some attributes stored in the simulation state become anonymous attributes in geometry nodes.
+   * This maps attribute names to their corresponding anonymous attribute ids. */
+  Map<std::string, AnonymousAttributeIDPtr> attribute_map;
+  Vector<GeometrySet *> geometries;
+
+  for (const int i : node_simulation_items.index_range()) {
+    const NodeSimulationItem &item = node_simulation_items[i];
+    const eNodeSocketDatatype socket_type = eNodeSocketDatatype(item.socket_type);
+    const CPPType &cpp_type = get_simulation_item_cpp_type(socket_type);
+
+    void *r_output_value = r_output_values[i];
+
+    if (!zone_state.item_by_identifier.contains(item.identifier)) {
+      cpp_type.value_initialize(r_output_value);
+      continue;
+    }
+    const bke::sim::SimulationStateItem &state_item = *zone_state.item_by_identifier.lookup(
+        item.identifier);
+
+    switch (socket_type) {
+      case SOCK_GEOMETRY: {
+        if (const auto *geo_state_item =
+                dynamic_cast<const bke::sim::GeometrySimulationStateItem *>(&state_item))
+        {
+          GeometrySet *geometry = new (r_output_value) GeometrySet(geo_state_item->geometry());
+          geometries.append(geometry);
+        }
+        else {
+          cpp_type.value_initialize(r_output_value);
+        }
+        break;
+      }
+      case SOCK_FLOAT:
+      case SOCK_VECTOR:
+      case SOCK_INT:
+      case SOCK_BOOLEAN:
+      case SOCK_RGBA: {
+        const fn::ValueOrFieldCPPType &value_or_field_type =
+            *fn::ValueOrFieldCPPType::get_from_self(cpp_type);
+        if (const auto *primitive_state_item =
+                dynamic_cast<const bke::sim::PrimitiveSimulationStateItem *>(&state_item))
+        {
+          if (primitive_state_item->type() == value_or_field_type.value) {
+            value_or_field_type.construct_from_value(r_output_value,
+                                                     primitive_state_item->value());
+          }
+          else {
+            cpp_type.value_initialize(r_output_value);
+          }
+        }
+        else if (const auto *attribute_state_item =
+                     dynamic_cast<const bke::sim::AttributeSimulationStateItem *>(&state_item))
+        {
+          AnonymousAttributeIDPtr attribute_id = MEM_new<NodeAnonymousAttributeID>(
+              __func__,
+              self_object,
+              compute_context,
+              node,
+              std::to_string(item.identifier),
+              item.name);
+          GField field{std::make_shared<AnonymousAttributeFieldInput>(
+              attribute_id, value_or_field_type.value, node.label_or_name())};
+          value_or_field_type.construct_from_field(r_output_value, std::move(field));
+          attribute_map.add(attribute_state_item->name(), std::move(attribute_id));
+        }
+        else {
+          cpp_type.value_initialize(r_output_value);
+        }
+        break;
+      }
+      case SOCK_STRING: {
+        if (const auto *string_state_item =
+                dynamic_cast<const bke::sim::StringSimulationStateItem *>(&state_item))
+        {
+          new (r_output_value) ValueOrField<std::string>(string_state_item->value());
+        }
+        else {
+          cpp_type.value_initialize(r_output_value);
+        }
+        break;
+      }
+      default: {
+        cpp_type.value_initialize(r_output_value);
+        break;
+      }
+    }
+  }
+
+  /* Make some attributes anonymous. */
+  for (GeometrySet *geometry : geometries) {
+    for (const GeometryComponentType type : {GEO_COMPONENT_TYPE_MESH,
+                                             GEO_COMPONENT_TYPE_CURVE,
+                                             GEO_COMPONENT_TYPE_POINT_CLOUD,
+                                             GEO_COMPONENT_TYPE_INSTANCES})
+    {
+      if (!geometry->has(type)) {
+        continue;
+      }
+      GeometryComponent &component = geometry->get_component_for_write(type);
+      MutableAttributeAccessor attributes = *component.attributes_for_write();
+      for (const MapItem<std::string, AnonymousAttributeIDPtr> &attribute_item :
+           attribute_map.items()) {
+        attributes.rename(attribute_item.key, *attribute_item.value);
+      }
+    }
+  }
+}
+
+void values_to_simulation_state(const Span<NodeSimulationItem> node_simulation_items,
+                                const Span<void *> input_values,
+                                bke::sim::SimulationZoneState &r_zone_state)
+{
+  Vector<GeometrySet *> stored_geometries;
+
+  for (const int i : node_simulation_items.index_range()) {
+    const NodeSimulationItem &item = node_simulation_items[i];
+    const eNodeSocketDatatype socket_type = eNodeSocketDatatype(item.socket_type);
+    void *input_value = input_values[i];
+
+    std::unique_ptr<bke::sim::SimulationStateItem> state_item;
+    switch (socket_type) {
+      case SOCK_GEOMETRY: {
+        GeometrySet &geometry = *static_cast<GeometrySet *>(input_value);
+        auto geometry_state_item = std::make_unique<bke::sim::GeometrySimulationStateItem>(
+            std::move(geometry));
+        stored_geometries.append(&geometry_state_item->geometry());
+        state_item = std::move(geometry_state_item);
+        break;
+      }
+      case SOCK_FLOAT:
+      case SOCK_VECTOR:
+      case SOCK_INT:
+      case SOCK_BOOLEAN:
+      case SOCK_RGBA: {
+        const CPPType &type = get_simulation_item_cpp_type(item);
+        const fn::ValueOrFieldCPPType &value_or_field_type =
+            *fn::ValueOrFieldCPPType::get_from_self(type);
+        if (value_or_field_type.is_field(input_value)) {
+          /* Fields are evaluated and stored as attributes. */
+          if (!stored_geometries.is_empty()) {
+            /* Possible things to consider:
+             * - Store attributes on multiple/all geometries.
+             * - If the attribute is an anonymous attribute, just rename it for the simulation
+             *   state, without considering the domain. This would allow e.g. having the attribute
+             *   only on some parts of the geometry set.
+             */
+            GeometrySet &geometry = *stored_geometries.last();
+            const GField &field = *value_or_field_type.get_field_ptr(input_value);
+            const eAttrDomain domain = eAttrDomain(item.attribute_domain);
+            const std::string attribute_name = ".sim_" + std::to_string(item.identifier);
+            if (geometry.has_pointcloud()) {
+              PointCloudComponent &component =
+                  geometry.get_component_for_write<PointCloudComponent>();
+              bke::try_capture_field_on_geometry(component, attribute_name, domain, field);
+            }
+            if (geometry.has_mesh()) {
+              MeshComponent &component = geometry.get_component_for_write<MeshComponent>();
+              bke::try_capture_field_on_geometry(component, attribute_name, domain, field);
+            }
+            if (geometry.has_curves()) {
+              CurveComponent &component = geometry.get_component_for_write<CurveComponent>();
+              bke::try_capture_field_on_geometry(component, attribute_name, domain, field);
+            }
+            if (geometry.has_instances()) {
+              InstancesComponent &component =
+                  geometry.get_component_for_write<InstancesComponent>();
+              bke::try_capture_field_on_geometry(component, attribute_name, domain, field);
+            }
+            state_item = std::make_unique<bke::sim::AttributeSimulationStateItem>(attribute_name);
+          }
+        }
+        else {
+          const void *value = value_or_field_type.get_value_ptr(input_value);
+          state_item = std::make_unique<bke::sim::PrimitiveSimulationStateItem>(
+              value_or_field_type.value, value);
+        }
+        break;
+      }
+      case SOCK_STRING: {
+        const ValueOrField<std::string> &value = *static_cast<const ValueOrField<std::string> *>(
+            input_value);
+        state_item = std::make_unique<bke::sim::StringSimulationStateItem>(value.as_value());
+        break;
+      }
+      default:
+        break;
+    }
+
+    if (state_item) {
+      r_zone_state.item_by_identifier.add_new(item.identifier, std::move(state_item));
+    }
+  }
+
+  for (GeometrySet *geometry : stored_geometries) {
+    cleanup_geometry_for_simulation_state(*geometry);
+    geometry->ensure_owns_all_data();
+  }
+}
+
+}  // namespace blender::nodes
+
+namespace blender::nodes::node_geo_simulation_output_cc {
+
+NODE_STORAGE_FUNCS(NodeGeometrySimulationOutput);
+
+struct EvalData {
+  bool is_first_evaluation = true;
+};
+
+class LazyFunctionForSimulationOutputNode final : public LazyFunction {
+  const bNode &node_;
+  Span<NodeSimulationItem> simulation_items_;
+
+ public:
+  LazyFunctionForSimulationOutputNode(const bNode &node,
+                                      GeometryNodesLazyFunctionGraphInfo &own_lf_graph_info)
+      : node_(node)
+  {
+    debug_name_ = "Simulation Output";
+    const NodeGeometrySimulationOutput &storage = node_storage(node);
+    simulation_items_ = {storage.items, storage.items_num};
+
+    MutableSpan<int> lf_index_by_bsocket = own_lf_graph_info.mapping.lf_index_by_bsocket;
+
+    for (const int i : simulation_items_.index_range()) {
+      const NodeSimulationItem &item = simulation_items_[i];
+      const bNodeSocket &input_bsocket = node.input_socket(i);
+      const bNodeSocket &output_bsocket = node.output_socket(i);
+
+      const CPPType &type = get_simulation_item_cpp_type(item);
+
+      lf_index_by_bsocket[input_bsocket.index_in_tree()] = inputs_.append_and_get_index_as(
+          item.name, type, lf::ValueUsage::Maybe);
+      lf_index_by_bsocket[output_bsocket.index_in_tree()] = outputs_.append_and_get_index_as(
+          item.name, type);
+    }
+  }
+
+  void *init_storage(LinearAllocator<> &allocator) const
+  {
+    return allocator.construct<EvalData>().get();
+  }
+
+  void destruct_storage(void *storage) const
+  {
+    std::destroy_at(static_cast<EvalData *>(storage));
+  }
+
+  void execute_impl(lf::Params &params, const lf::Context &context) const final
+  {
+    GeoNodesLFUserData &user_data = *static_cast<GeoNodesLFUserData *>(context.user_data);
+    GeoNodesModifierData &modifier_data = *user_data.modifier_data;
+    EvalData &eval_data = *static_cast<EvalData *>(context.storage);
+    BLI_SCOPED_DEFER([&]() { eval_data.is_first_evaluation = false; });
+
+    const bke::sim::SimulationZoneID zone_id = get_simulation_zone_id(*user_data.compute_context,
+                                                                      node_.identifier);
+
+    const bke::sim::SimulationZoneState *current_zone_state =
+        modifier_data.current_simulation_state ?
+            modifier_data.current_simulation_state->get_zone_state(zone_id) :
+            nullptr;
+    if (eval_data.is_first_evaluation && current_zone_state != nullptr) {
+      /* Common case when data is cached already. */
+      this->output_cached_state(
+          params, *modifier_data.self_object, *user_data.compute_context, *current_zone_state);
+      return;
+    }
+
+    if (modifier_data.current_simulation_state_for_write == nullptr) {
+      const bke::sim::SimulationZoneState *prev_zone_state =
+          modifier_data.prev_simulation_state ?
+              modifier_data.prev_simulation_state->get_zone_state(zone_id) :
+              nullptr;
+      if (prev_zone_state == nullptr) {
+        /* There is no previous simulation state and we also don't create a new one, so just output
+         * defaults. */
+        params.set_default_remaining_outputs();
+        return;
+      }
+      const bke::sim::SimulationZoneState *next_zone_state =
+          modifier_data.next_simulation_state ?
+              modifier_data.next_simulation_state->get_zone_state(zone_id) :
+              nullptr;
+      if (next_zone_state == nullptr) {
+        /* Output the last cached simulation state. */
+        this->output_cached_state(
+            params, *modifier_data.self_object, *user_data.compute_context, *prev_zone_state);
+        return;
+      }
+      /* A previous and next frame is cached already, but the current frame is not. */
+      this->output_mixed_cached_state(params,
+                                      *modifier_data.self_object,
+                                      *user_data.compute_context,
+                                      *prev_zone_state,
+                                      *next_zone_state,
+                                      modifier_data.simulation_state_mix_factor);
+      return;
+    }
+
+    bke::sim::SimulationZoneState &new_zone_state =
+        modifier_data.current_simulation_state_for_write->get_zone_state_for_write(zone_id);
+    if (eval_data.is_first_evaluation) {
+      new_zone_state.item_by_identifier.clear();
+    }
+
+    Array<void *> input_values(simulation_items_.size(), nullptr);
+    for (const int i : simulation_items_.index_range()) {
+      input_values[i] = params.try_get_input_data_ptr_or_request(i);
+    }
+    if (input_values.as_span().contains(nullptr)) {
+      /* Wait until all inputs are available. */
+      return;
+    }
+    values_to_simulation_state(simulation_items_, input_values, new_zone_state);
+    this->output_cached_state(
+        params, *modifier_data.self_object, *user_data.compute_context, new_zone_state);
+  }
+
+  void output_cached_state(lf::Params &params,
+                           const Object &self_object,
+                           const ComputeContext &compute_context,
+                           const bke::sim::SimulationZoneState &state) const
+  {
+    Array<void *> output_values(simulation_items_.size());
+    for (const int i : simulation_items_.index_range()) {
+      output_values[i] = params.get_output_data_ptr(i);
+    }
+    simulation_state_to_values(
+        simulation_items_, state, self_object, compute_context, node_, output_values);
+    for (const int i : simulation_items_.index_range()) {
+      params.output_set(i);
+    }
+  }
+
+  void output_mixed_cached_state(lf::Params &params,
+                                 const Object &self_object,
+                                 const ComputeContext &compute_context,
+                                 const bke::sim::SimulationZoneState &prev_state,
+                                 const bke::sim::SimulationZoneState &next_state,
+                                 const float mix_factor) const
+  {
+    /* TODO: Implement subframe mixing. */
+    this->output_cached_state(params, self_object, compute_context, prev_state);
+    UNUSED_VARS(next_state, mix_factor);
+  }
+};
+
+}  // namespace blender::nodes::node_geo_simulation_output_cc
+
+namespace blender::nodes {
+
+std::unique_ptr<LazyFunction> get_simulation_output_lazy_function(
+    const bNode &node, GeometryNodesLazyFunctionGraphInfo &own_lf_graph_info)
+{
+  namespace file_ns = blender::nodes::node_geo_simulation_output_cc;
+  BLI_assert(node.type == GEO_NODE_SIMULATION_OUTPUT);
+  return std::make_unique<file_ns::LazyFunctionForSimulationOutputNode>(node, own_lf_graph_info);
+}
+
+bke::sim::SimulationZoneID get_simulation_zone_id(const ComputeContext &compute_context,
+                                                  const int output_node_id)
+{
+  bke::sim::SimulationZoneID zone_id;
+  for (const ComputeContext *context = &compute_context; context != nullptr;
+       context = context->parent())
+  {
+    if (const auto *node_context = dynamic_cast<const bke::NodeGroupComputeContext *>(context)) {
+      zone_id.node_ids.append(node_context->node_id());
+    }
+  }
+  std::reverse(zone_id.node_ids.begin(), zone_id.node_ids.end());
+  zone_id.node_ids.append(output_node_id);
+  return zone_id;
+}
+
+}  // namespace blender::nodes
+
+namespace blender::nodes::node_geo_simulation_output_cc {
+
+static void node_declare_dynamic(const bNodeTree & /*node_tree*/,
+                                 const bNode &node,
+                                 NodeDeclaration &r_declaration)
+{
+  const NodeGeometrySimulationOutput &storage = node_storage(node);
+  socket_declarations_for_simulation_items({storage.items, storage.items_num}, r_declaration);
+}
+
+static void search_node_add_ops(GatherAddNodeSearchParams &params)
+{
+  AddNodeItem item;
+  item.ui_name = IFACE_("Simulation Zone");
+  item.description = TIP_("Add a new simulation input and output nodes to the node tree");
+  item.add_fn = [](const bContext &C, bNodeTree &node_tree, float2 cursor) {
+    bNode *input = nodeAddNode(&C, &node_tree, "GeometryNodeSimulationInput");
+    bNode *output = nodeAddNode(&C, &node_tree, "GeometryNodeSimulationOutput");
+    static_cast<NodeGeometrySimulationInput *>(input->storage)->output_node_id =
+        output->identifier;
+
+    NodeSimulationItem &item = node_storage(*output).items[0];
+
+    update_node_declaration_and_sockets(node_tree, *input);
+    update_node_declaration_and_sockets(node_tree, *output);
+
+    nodeAddLink(
+        &node_tree,
+        input,
+        nodeFindSocket(input, SOCK_OUT, socket_identifier_for_simulation_item(item).c_str()),
+        output,
+        nodeFindSocket(output, SOCK_IN, socket_identifier_for_simulation_item(item).c_str()));
+
+    input->locx = cursor.x / UI_SCALE_FAC - 150;
+    input->locy = cursor.y / UI_SCALE_FAC + 20;
+    output->locx = cursor.x / UI_SCALE_FAC + 150;
+    output->locy = cursor.y / UI_SCALE_FAC + 20;
+
+    return Vector<bNode *>({input, output});
+  };
+  params.add_item(std::move(item));
+}
+
+static void node_init(bNodeTree * /*tree*/, bNode *node)
+{
+  NodeGeometrySimulationOutput *data = MEM_cnew<NodeGeometrySimulationOutput>(__func__);
+
+  data->next_identifier = 0;
+
+  data->items = MEM_cnew_array<NodeSimulationItem>(1, __func__);
+  data->items[0].name = BLI_strdup(DATA_("Geometry"));
+  data->items[0].socket_type = SOCK_GEOMETRY;
+  data->items[0].identifier = data->next_identifier++;
+  data->items_num = 1;
+
+  node->storage = data;
+}
+
+static void node_free_storage(bNode *node)
+{
+  if (!node->storage) {
+    return;
+  }
+  NodeGeometrySimulationOutput &storage = node_storage(*node);
+  for (NodeSimulationItem &item : MutableSpan(storage.items, storage.items_num)) {
+    MEM_SAFE_FREE(item.name);
+  }
+  MEM_SAFE_FREE(storage.items);
+  MEM_freeN(node->storage);
+}
+
+static void node_copy_storage(bNodeTree * /*dst_tree*/, bNode *dst_node, const bNode *src_node)
+{
+  const NodeGeometrySimulationOutput &src_storage = node_storage(*src_node);
+  NodeGeometrySimulationOutput *dst_storage = MEM_cnew<NodeGeometrySimulationOutput>(__func__);
+
+  dst_storage->items = MEM_cnew_array<NodeSimulationItem>(src_storage.items_num, __func__);
+  dst_storage->items_num = src_storage.items_num;
+  dst_storage->active_index = src_storage.active_index;
+  dst_storage->next_identifier = src_storage.next_identifier;
+  for (const int i : IndexRange(src_storage.items_num)) {
+    if (char *name = src_storage.items[i].name) {
+      dst_storage->items[i].identifier = src_storage.items[i].identifier;
+      dst_storage->items[i].name = BLI_strdup(name);
+      dst_storage->items[i].socket_type = src_storage.items[i].socket_type;
+      dst_storage->items[i].attribute_domain = src_storage.items[i].attribute_domain;
+    }
+  }
+
+  dst_node->storage = dst_storage;
+}
+
+static bool node_insert_link(bNodeTree *ntree, bNode *node, bNodeLink *link)
+{
+  NodeGeometrySimulationOutput &storage = node_storage(*node);
+  if (link->tonode == node) {
+    if (link->tosock->identifier == StringRef("__extend__")) {
+      if (const NodeSimulationItem *item = NOD_geometry_simulation_output_add_item_from_socket(
+              &storage, link->fromnode, link->fromsock))
+      {
+        update_node_declaration_and_sockets(*ntree, *node);
+        link->tosock = nodeFindSocket(
+            node, SOCK_IN, socket_identifier_for_simulation_item(*item).c_str());
+      }
+      else {
+        return false;
+      }
+    }
+  }
+  else {
+    BLI_assert(link->fromnode == node);
+    if (link->fromsock->identifier == StringRef("__extend__")) {
+      if (const NodeSimulationItem *item = NOD_geometry_simulation_output_add_item_from_socket(
+              &storage, link->fromnode, link->tosock))
+      {
+        update_node_declaration_and_sockets(*ntree, *node);
+        link->fromsock = nodeFindSocket(
+            node, SOCK_OUT, socket_identifier_for_simulation_item(*item).c_str());
+      }
+      else {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+}  // namespace blender::nodes::node_geo_simulation_output_cc
+
+void register_node_type_geo_simulation_output()
+{
+  namespace file_ns = blender::nodes::node_geo_simulation_output_cc;
+
+  static bNodeType ntype;
+
+  geo_node_type_base(
+      &ntype, GEO_NODE_SIMULATION_OUTPUT, "Simulation Output", NODE_CLASS_INTERFACE);
+  ntype.initfunc = file_ns::node_init;
+  ntype.declare_dynamic = file_ns::node_declare_dynamic;
+  ntype.gather_add_node_search_ops = file_ns::search_node_add_ops;
+  ntype.gather_link_search_ops = nullptr;
+  ntype.insert_link = file_ns::node_insert_link;
+  node_type_storage(&ntype,
+                    "NodeGeometrySimulationOutput",
+                    file_ns::node_free_storage,
+                    file_ns::node_copy_storage);
+  nodeRegisterType(&ntype);
+}
+
+blender::Span<NodeSimulationItem> NodeGeometrySimulationOutput::items_span() const
+{
+  return blender::Span<NodeSimulationItem>(items, items_num);
+}
+
+blender::MutableSpan<NodeSimulationItem> NodeGeometrySimulationOutput::items_span_for_write()
+{
+  return blender::MutableSpan<NodeSimulationItem>(items, items_num);
+}
+
+blender::IndexRange NodeGeometrySimulationOutput::items_range() const
+{
+  return blender::IndexRange(items_num);
+}
+
+bool NOD_geometry_simulation_output_item_socket_type_supported(
+    const eNodeSocketDatatype socket_type)
+{
+  return ELEM(socket_type,
+              SOCK_FLOAT,
+              SOCK_VECTOR,
+              SOCK_RGBA,
+              SOCK_BOOLEAN,
+              SOCK_INT,
+              SOCK_STRING,
+              SOCK_GEOMETRY);
+}
+
+bNode *NOD_geometry_simulation_output_find_node_by_item(bNodeTree *ntree,
+                                                        const NodeSimulationItem *item)
+{
+  ntree->ensure_topology_cache();
+  for (bNode *node : ntree->nodes_by_type("GeometryNodeSimulationOutput")) {
+    NodeGeometrySimulationOutput *sim = static_cast<NodeGeometrySimulationOutput *>(node->storage);
+    if (sim->items_span().contains_ptr(item)) {
+      return node;
+    }
+  }
+  return nullptr;
+}
+
+bool NOD_geometry_simulation_output_item_set_unique_name(NodeGeometrySimulationOutput *sim,
+                                                         NodeSimulationItem *item,
+                                                         const char *name,
+                                                         const char *defname)
+{
+  char unique_name[MAX_NAME + 4];
+  BLI_strncpy(unique_name, name, sizeof(unique_name));
+
+  blender::nodes::SimulationItemsUniqueNameArgs args{sim, item};
+  const bool name_changed = BLI_uniquename_cb(blender::nodes::simulation_items_unique_name_check,
+                                              &args,
+                                              defname,
+                                              '.',
+                                              unique_name,
+                                              ARRAY_SIZE(unique_name));
+  item->name = BLI_strdup(unique_name);
+  return name_changed;
+}
+
+bool NOD_geometry_simulation_output_contains_item(NodeGeometrySimulationOutput *sim,
+                                                  const NodeSimulationItem *item)
+{
+  return sim->items_span().contains_ptr(item);
+}
+
+NodeSimulationItem *NOD_geometry_simulation_output_get_active_item(
+    NodeGeometrySimulationOutput *sim)
+{
+  if (!sim->items_range().contains(sim->active_index)) {
+    return nullptr;
+  }
+  return &sim->items[sim->active_index];
+}
+
+void NOD_geometry_simulation_output_set_active_item(NodeGeometrySimulationOutput *sim,
+                                                    NodeSimulationItem *item)
+{
+  if (sim->items_span().contains_ptr(item)) {
+    sim->active_index = item - sim->items;
+  }
+}
+
+NodeSimulationItem *NOD_geometry_simulation_output_find_item(NodeGeometrySimulationOutput *sim,
+                                                             const char *name)
+{
+  for (NodeSimulationItem &item : sim->items_span_for_write()) {
+    if (STREQ(item.name, name)) {
+      return &item;
+    }
+  }
+  return nullptr;
+}
+
+NodeSimulationItem *NOD_geometry_simulation_output_add_item(NodeGeometrySimulationOutput *sim,
+                                                            const short socket_type,
+                                                            const char *name)
+{
+  return NOD_geometry_simulation_output_insert_item(sim, socket_type, name, sim->items_num);
+}
+
+NodeSimulationItem *NOD_geometry_simulation_output_insert_item(NodeGeometrySimulationOutput *sim,
+                                                               const short socket_type,
+                                                               const char *name,
+                                                               int index)
+{
+  if (!NOD_geometry_simulation_output_item_socket_type_supported(eNodeSocketDatatype(socket_type)))
+  {
+    return nullptr;
+  }
+
+  NodeSimulationItem *old_items = sim->items;
+  sim->items = MEM_cnew_array<NodeSimulationItem>(sim->items_num + 1, __func__);
+  for (const int i : blender::IndexRange(index)) {
+    sim->items[i] = old_items[i];
+  }
+  for (const int i : blender::IndexRange(index, sim->items_num - index)) {
+    sim->items[i + 1] = old_items[i];
+  }
+
+  const char *defname = nodeStaticSocketLabel(socket_type, 0);
+  NodeSimulationItem &added_item = sim->items[index];
+  added_item.identifier = sim->next_identifier++;
+  NOD_geometry_simulation_output_item_set_unique_name(sim, &added_item, name, defname);
+  added_item.socket_type = socket_type;
+
+  sim->items_num++;
+  MEM_SAFE_FREE(old_items);
+
+  return &added_item;
+}
+
+NodeSimulationItem *NOD_geometry_simulation_output_add_item_from_socket(
+    NodeGeometrySimulationOutput *sim, const bNode * /*from_node*/, const bNodeSocket *from_sock)
+{
+  return NOD_geometry_simulation_output_insert_item(
+      sim, from_sock->type, from_sock->name, sim->items_num);
+}
+
+NodeSimulationItem *NOD_geometry_simulation_output_insert_item_from_socket(
+    NodeGeometrySimulationOutput *sim,
+    const bNode * /*from_node*/,
+    const bNodeSocket *from_sock,
+    int index)
+{
+  return NOD_geometry_simulation_output_insert_item(sim, from_sock->type, from_sock->name, index);
+}
+
+void NOD_geometry_simulation_output_remove_item(NodeGeometrySimulationOutput *sim,
+                                                NodeSimulationItem *item)
+{
+  const int index = item - sim->items;
+  if (index < 0 || index >= sim->items_num) {
+    return;
+  }
+
+  NodeSimulationItem *old_items = sim->items;
+  sim->items = MEM_cnew_array<NodeSimulationItem>(sim->items_num - 1, __func__);
+  for (const int i : blender::IndexRange(index)) {
+    sim->items[i] = old_items[i];
+  }
+  for (const int i : blender::IndexRange(index, sim->items_num - index).drop_front(1)) {
+    sim->items[i - 1] = old_items[i];
+  }
+
+  MEM_SAFE_FREE(old_items[index].name);
+
+  sim->items_num--;
+  MEM_SAFE_FREE(old_items);
+}
+
+void NOD_geometry_simulation_output_clear_items(struct NodeGeometrySimulationOutput *sim)
+{
+  for (NodeSimulationItem &item : sim->items_span_for_write()) {
+    MEM_SAFE_FREE(item.name);
+  }
+  MEM_SAFE_FREE(sim->items);
+  sim->items = nullptr;
+  sim->items_num = 0;
+}
+
+void NOD_geometry_simulation_output_move_item(NodeGeometrySimulationOutput *sim,
+                                              int from_index,
+                                              int to_index)
+{
+  BLI_assert(from_index >= 0 && from_index < sim->items_num);
+  BLI_assert(to_index >= 0 && to_index < sim->items_num);
+
+  if (from_index == to_index) {
+    return;
+  }
+
+  if (from_index < to_index) {
+    const NodeSimulationItem tmp = sim->items[from_index];
+    for (int i = from_index; i < to_index; ++i) {
+      sim->items[i] = sim->items[i + 1];
+    }
+    sim->items[to_index] = tmp;
+  }
+  else /* from_index > to_index */ {
+    const NodeSimulationItem tmp = sim->items[from_index];
+    for (int i = from_index; i > to_index; --i) {
+      sim->items[i] = sim->items[i - 1];
+    }
+    sim->items[to_index] = tmp;
+  }
+}

source/blender/nodes/intern/geometry_nodes_lazy_function.cc

@@ -105,42 +105,33 @@ static void lazy_function_interface_from_node(const bNode &node,
   }
 }
 
-/**
- * An anonymous attribute created by a node.
- */
-class NodeAnonymousAttributeID : public AnonymousAttributeID {
-  std::string long_name_;
-  std::string socket_name_;
-
- public:
-  NodeAnonymousAttributeID(const Object &object,
-                           const ComputeContext &compute_context,
-                           const bNode &bnode,
-                           const StringRef identifier,
-                           const StringRef name)
-      : socket_name_(name)
+NodeAnonymousAttributeID::NodeAnonymousAttributeID(const Object &object,
+                                                   const ComputeContext &compute_context,
+                                                   const bNode &bnode,
+                                                   const StringRef identifier,
+                                                   const StringRef name)
+    : socket_name_(name)
+{
+  const ComputeContextHash &hash = compute_context.hash();
   {
-    const ComputeContextHash &hash = compute_context.hash();
-    {
-      std::stringstream ss;
-      ss << hash << "_" << object.id.name << "_" << bnode.identifier << "_" << identifier;
-      long_name_ = ss.str();
-    }
-    {
-      uint64_t hash_result[2];
-      BLI_hash_md5_buffer(long_name_.data(), long_name_.size(), hash_result);
-      std::stringstream ss;
-      ss << ".a_" << std::hex << hash_result[0] << hash_result[1];
-      name_ = ss.str();
-      BLI_assert(name_.size() < MAX_CUSTOMDATA_LAYER_NAME);
-    }
+    std::stringstream ss;
+    ss << hash << "_" << object.id.name << "_" << bnode.identifier << "_" << identifier;
+    long_name_ = ss.str();
   }
-
-  std::string user_name() const override
   {
-    return socket_name_;
+    uint64_t hash_result[2];
+    BLI_hash_md5_buffer(long_name_.data(), long_name_.size(), hash_result);
+    std::stringstream ss;
+    ss << ".a_" << std::hex << hash_result[0] << hash_result[1];
+    name_ = ss.str();
+    BLI_assert(name_.size() < MAX_CUSTOMDATA_LAYER_NAME);
   }
-};
+}
+
+std::string NodeAnonymousAttributeID::user_name() const
+{
+  return socket_name_;
+}
 
 /**
  * Used for most normal geometry nodes like Subdivision Surface and Set Position.
@@ -846,6 +837,31 @@ class LazyFunctionForViewerInputUsage : public LazyFunction {
   }
 };
 
+class LazyFunctionForSimulationInputsUsage : public LazyFunction {
+ private:
+  const bNode &sim_output_node_;
+
+ public:
+  LazyFunctionForSimulationInputsUsage(const bNode &sim_output_node)
+      : sim_output_node_(sim_output_node)
+  {
+    debug_name_ = "Simulation Inputs Usage";
+    outputs_.append_as("Is Initialization", CPPType::get<bool>());
+    outputs_.append_as("Do Simulation Step", CPPType::get<bool>());
+  }
+
+  void execute_impl(lf::Params &params, const lf::Context &context) const override
+  {
+    const GeoNodesLFUserData &user_data = *static_cast<GeoNodesLFUserData *>(context.user_data);
+    const GeoNodesModifierData &modifier_data = *user_data.modifier_data;
+
+    params.set_output(0,
+                      modifier_data.current_simulation_state_for_write != nullptr &&
+                          modifier_data.prev_simulation_state == nullptr);
+    params.set_output(1, modifier_data.current_simulation_state_for_write != nullptr);
+  }
+};
+
 /**
  * This lazy-function wraps a group node. Internally it just executes the lazy-function graph of
  * the referenced group.
@@ -1370,6 +1386,10 @@ struct GeometryNodesLazyFunctionGraphBuilder {
    * All group input nodes are combined into one dummy node in the lazy-function graph.
    */
   lf::DummyNode *group_input_lf_node_;
+  /**
+   * A #LazyFunctionForSimulationInputsUsage for each simulation zone.
+   */
+  Map<const bNode *, lf::Node *> simulation_inputs_usage_nodes_;
 
   friend class UsedSocketVisualizeOptions;
 
@@ -1508,6 +1528,14 @@ struct GeometryNodesLazyFunctionGraphBuilder {
           this->handle_viewer_node(*bnode);
           break;
         }
+        case GEO_NODE_SIMULATION_INPUT: {
+          this->handle_simulation_input_node(btree_, *bnode);
+          break;
+        }
+        case GEO_NODE_SIMULATION_OUTPUT: {
+          this->handle_simulation_output_node(*bnode);
+          break;
+        }
         case GEO_NODE_SWITCH: {
           this->handle_switch_node(*bnode);
           break;
@@ -1785,6 +1813,60 @@ struct GeometryNodesLazyFunctionGraphBuilder {
     mapping_->viewer_node_map.add(&bnode, &lf_node);
   }
 
+  void handle_simulation_input_node(const bNodeTree &node_tree, const bNode &bnode)
+  {
+    const NodeGeometrySimulationInput *storage = static_cast<const NodeGeometrySimulationInput *>(
+        bnode.storage);
+    if (node_tree.node_by_id(storage->output_node_id) == nullptr) {
+      return;
+    }
+
+    std::unique_ptr<LazyFunction> lazy_function = get_simulation_input_lazy_function(
+        node_tree, bnode, *lf_graph_info_);
+    lf::FunctionNode &lf_node = lf_graph_->add_function(*lazy_function);
+    lf_graph_info_->functions.append(std::move(lazy_function));
+
+    for (const int i : bnode.input_sockets().index_range().drop_back(1)) {
+      const bNodeSocket &bsocket = bnode.input_socket(i);
+      lf::InputSocket &lf_socket = lf_node.input(
+          mapping_->lf_index_by_bsocket[bsocket.index_in_tree()]);
+      input_socket_map_.add(&bsocket, &lf_socket);
+      mapping_->bsockets_by_lf_socket_map.add(&lf_socket, &bsocket);
+    }
+    for (const int i : bnode.output_sockets().index_range().drop_back(1)) {
+      const bNodeSocket &bsocket = bnode.output_socket(i);
+      lf::OutputSocket &lf_socket = lf_node.output(
+          mapping_->lf_index_by_bsocket[bsocket.index_in_tree()]);
+      output_socket_map_.add(&bsocket, &lf_socket);
+      mapping_->bsockets_by_lf_socket_map.add(&lf_socket, &bsocket);
+    }
+  }
+
+  void handle_simulation_output_node(const bNode &bnode)
+  {
+    std::unique_ptr<LazyFunction> lazy_function = get_simulation_output_lazy_function(
+        bnode, *lf_graph_info_);
+    lf::FunctionNode &lf_node = lf_graph_->add_function(*lazy_function);
+    lf_graph_info_->functions.append(std::move(lazy_function));
+
+    for (const int i : bnode.input_sockets().index_range().drop_back(1)) {
+      const bNodeSocket &bsocket = bnode.input_socket(i);
+      lf::InputSocket &lf_socket = lf_node.input(
+          mapping_->lf_index_by_bsocket[bsocket.index_in_tree()]);
+      input_socket_map_.add(&bsocket, &lf_socket);
+      mapping_->bsockets_by_lf_socket_map.add(&lf_socket, &bsocket);
+    }
+    for (const int i : bnode.output_sockets().index_range().drop_back(1)) {
+      const bNodeSocket &bsocket = bnode.output_socket(i);
+      lf::OutputSocket &lf_socket = lf_node.output(
+          mapping_->lf_index_by_bsocket[bsocket.index_in_tree()]);
+      output_socket_map_.add(&bsocket, &lf_socket);
+      mapping_->bsockets_by_lf_socket_map.add(&lf_socket, &bsocket);
+    }
+
+    mapping_->sim_output_node_map.add(&bnode, &lf_node);
+  }
+
   void handle_switch_node(const bNode &bnode)
   {
     std::unique_ptr<LazyFunction> lazy_function = get_switch_node_lazy_function(bnode);
@@ -2136,6 +2218,14 @@ struct GeometryNodesLazyFunctionGraphBuilder {
           this->build_viewer_node_socket_usage(*bnode);
           break;
         }
+        case GEO_NODE_SIMULATION_INPUT: {
+          this->build_simulation_input_socket_usage(*bnode);
+          break;
+        }
+        case GEO_NODE_SIMULATION_OUTPUT: {
+          this->build_simulation_output_socket_usage(*bnode);
+          break;
+        }
         case NODE_GROUP:
         case NODE_CUSTOM_GROUP: {
           this->build_group_node_socket_usage(*bnode, or_socket_usages_cache);
@@ -2296,6 +2386,44 @@ struct GeometryNodesLazyFunctionGraphBuilder {
     }
   }
 
+  void build_simulation_input_socket_usage(const bNode &bnode)
+  {
+    const NodeGeometrySimulationInput *storage = static_cast<const NodeGeometrySimulationInput *>(
+        bnode.storage);
+    const bNode *sim_output_node = btree_.node_by_id(storage->output_node_id);
+    if (sim_output_node == nullptr) {
+      return;
+    }
+    lf::Node &lf_node = this->get_simulation_inputs_usage_node(*sim_output_node);
+    for (const bNodeSocket *bsocket : bnode.input_sockets()) {
+      if (bsocket->is_available()) {
+        socket_is_used_map_[bsocket->index_in_tree()] = &lf_node.output(0);
+      }
+    }
+  }
+
+  void build_simulation_output_socket_usage(const bNode &bnode)
+  {
+    lf::Node &lf_node = this->get_simulation_inputs_usage_node(bnode);
+    for (const bNodeSocket *bsocket : bnode.input_sockets()) {
+      if (bsocket->is_available()) {
+        socket_is_used_map_[bsocket->index_in_tree()] = &lf_node.output(1);
+      }
+    }
+  }
+
+  lf::Node &get_simulation_inputs_usage_node(const bNode &sim_output_bnode)
+  {
+    BLI_assert(sim_output_bnode.type == GEO_NODE_SIMULATION_OUTPUT);
+    return *simulation_inputs_usage_nodes_.lookup_or_add_cb(&sim_output_bnode, [&]() {
+      auto lazy_function = std::make_unique<LazyFunctionForSimulationInputsUsage>(
+          sim_output_bnode);
+      lf::Node &lf_node = lf_graph_->add_function(*lazy_function);
+      lf_graph_info_->functions.append(std::move(lazy_function));
+      return &lf_node;
+    });
+  }
+
   void build_group_node_socket_usage(const bNode &bnode,
                                      OrSocketUsagesCache &or_socket_usages_cache)
   {

source/blender/windowmanager/WM_api.h

@@ -1516,6 +1516,7 @@ typedef enum eWM_JobType {
   WM_JOB_TYPE_LINEART,
   WM_JOB_TYPE_SEQ_DRAW_THUMBNAIL,
   WM_JOB_TYPE_SEQ_DRAG_DROP_PREVIEW,
+  WM_JOB_TYPE_BAKE_SIMULATION_NODES,
   /* add as needed, bake, seq proxy build
    * if having hard coded values is a problem */
 } eWM_JobType;