Walls, Flow, and Microfluidics 🧬

.pre-commit-config.yaml

@@ -0,0 +1,46 @@
+# Checks that run on every commit. They cover the fast, purely local gates:
+# formatting and linting for each language in the tree, plus the Python type
+# check. Tests, native builds, and backend conformance are not here; they need
+# a configured build and hardware, and run through CTest and the conformance
+# scripts instead.
+repos:
+  - repo: https://github.com/pre-commit/pre-commit-hooks
+    rev: v5.0.0
+    hooks:
+      - id: check-yaml
+      - id: check-toml
+      - id: check-merge-conflict
+      - id: check-added-large-files
+        # The device photomasks under docs/tutorials/devices are legitimately
+        # large and are tracked deliberately.
+        args: [--maxkb=1024]
+        exclude: ^docs/tutorials/devices/
+      - id: end-of-file-fixer
+        exclude: ^docs/tutorials/devices/
+      - id: trailing-whitespace
+        exclude: ^docs/tutorials/devices/
+
+  # Pinned to the version the project depends on, so a hook run and a local
+  # `ruff check` enforce the same rule set. Formatting is deliberately absent:
+  # the linter's rules, including import order, are what this tree follows.
+  - repo: https://github.com/astral-sh/ruff-pre-commit
+    rev: v0.16.3
+    hooks:
+      - id: ruff
+        args: [--fix]
+
+  - repo: https://github.com/pre-commit/mirrors-clang-format
+    rev: v19.1.7
+    hooks:
+      - id: clang-format
+        types_or: [c++, c, cuda, objective-c++, metal]
+        files: ^cpp/
+
+  - repo: local
+    hooks:
+      - id: pyright
+        name: pyright
+        entry: uv run pyright
+        language: system
+        pass_filenames: false
+        types: [python]

CMakeLists.txt

@@ -83,11 +83,12 @@ if(CM_ENABLE_METAL)
   add_custom_command(
     OUTPUT "${CM_METAL_COUPLED_RATES_HEADER}"
     COMMAND "${CMAKE_COMMAND}"
-      "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/coupled_rates.metal"
+      "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/grid_transport.metal;${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/coupled_rates.metal"
       "-DOUTPUT=${CM_METAL_COUPLED_RATES_HEADER}"
       "-DSYMBOL=coupled_rates_source"
       -P "${CMAKE_CURRENT_SOURCE_DIR}/cmake/EmbedMetalSource.cmake"
     DEPENDS
+      cpp/metal/kernels/grid_transport.metal
       cpp/metal/kernels/coupled_rates.metal
       cmake/EmbedMetalSource.cmake
     VERBATIM
@@ -107,11 +108,12 @@ if(CM_ENABLE_METAL)
   add_custom_command(
     OUTPUT "${CM_METAL_SIGNALS_HEADER}"
     COMMAND "${CMAKE_COMMAND}"
-      "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/signals.metal"
+      "-DINPUT=${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/grid_transport.metal;${CMAKE_CURRENT_SOURCE_DIR}/cpp/metal/kernels/signals.metal"
       "-DOUTPUT=${CM_METAL_SIGNALS_HEADER}"
       "-DSYMBOL=signals_source"
       -P "${CMAKE_CURRENT_SOURCE_DIR}/cmake/EmbedMetalSource.cmake"
     DEPENDS
+      cpp/metal/kernels/grid_transport.metal
       cpp/metal/kernels/signals.metal
       cmake/EmbedMetalSource.cmake
     VERBATIM

cmake/EmbedMetalSource.cmake

@@ -1,15 +1,23 @@
+# INPUT names one Metal source, or several to concatenate in order. A Metal
+# library is compiled from source at runtime with no include path, so a source
+# that shares helpers with another receives them by concatenation here.
 if(NOT DEFINED INPUT OR NOT DEFINED OUTPUT OR NOT DEFINED SYMBOL)
   message(FATAL_ERROR "EmbedMetalSource.cmake requires INPUT, OUTPUT, and SYMBOL")
 endif()
 if(NOT SYMBOL MATCHES "^[A-Za-z_][A-Za-z0-9_]*$")
   message(FATAL_ERROR "EmbedMetalSource.cmake received an invalid C++ symbol")
 endif()
 
-file(READ "${INPUT}" CM_METAL_SOURCE)
 get_filename_component(CM_OUTPUT_DIRECTORY "${OUTPUT}" DIRECTORY)
 file(MAKE_DIRECTORY "${CM_OUTPUT_DIRECTORY}")
 file(WRITE "${OUTPUT}" "#pragma once\n\nnamespace cm::metal {\ninline constexpr char ${SYMBOL}[] = R\"CM_METAL(")
-file(APPEND "${OUTPUT}" "${CM_METAL_SOURCE}")
+foreach(CM_METAL_INPUT IN LISTS INPUT)
+  file(READ "${CM_METAL_INPUT}" CM_METAL_SOURCE)
+  if(CM_METAL_SOURCE MATCHES "CM_METAL\\(" OR CM_METAL_SOURCE MATCHES "\\)CM_METAL")
+    message(FATAL_ERROR "Metal source ${CM_METAL_INPUT} contains the embedding delimiter")
+  endif()
+  file(APPEND "${OUTPUT}" "${CM_METAL_SOURCE}")
+endforeach()
 file(APPEND "${OUTPUT}" [=[)CM_METAL";
 }  // namespace cm::metal
 ]=])

cpp/core/constraints.cpp

@@ -3,6 +3,7 @@
 #include <cmath>
 #include <limits>
 #include <stdexcept>
+#include <string>
 #include <unordered_set>
 #include <utility>
 
@@ -36,45 +37,82 @@ void validate_sphere(const SphereConstraint& sphere) {
   }
   validate_coefficient(sphere.coefficient);
   switch (sphere.allowed_region) {
-    case SphereRegion::outside:
-    case SphereRegion::inside:
+    case ConstraintRegion::outside:
+    case ConstraintRegion::inside:
       return;
   }
   throw std::invalid_argument("checkpoint sphere uses an unknown allowed region");
 }
 
+bool positive_finite_extents(const Vec3& half_extents) {
+  return std::isfinite(half_extents.x) && half_extents.x > 0.0F && std::isfinite(half_extents.y) &&
+         half_extents.y > 0.0F && std::isfinite(half_extents.z) && half_extents.z > 0.0F;
+}
+
+void validate_box(const BoxConstraint& box) {
+  if (box.id == invalid_constraint_id || !finite(box.center) ||
+      !positive_finite_extents(box.half_extents)) {
+    throw std::invalid_argument("checkpoint box contains invalid geometry");
+  }
+  validate_coefficient(box.coefficient);
+  switch (box.allowed_region) {
+    case ConstraintRegion::outside:
+    case ConstraintRegion::inside:
+      return;
+  }
+  throw std::invalid_argument("checkpoint box uses an unknown allowed region");
+}
+
+void validate_cylinder(const CylinderConstraint& cylinder) {
+  if (cylinder.id == invalid_constraint_id || !finite(cylinder.center) ||
+      !std::isfinite(cylinder.radius) || cylinder.radius <= 0.0F ||
+      !std::isfinite(cylinder.half_height) || cylinder.half_height <= 0.0F) {
+    throw std::invalid_argument("checkpoint cylinder contains invalid geometry");
+  }
+  validate_coefficient(cylinder.coefficient);
+  switch (cylinder.allowed_region) {
+    case ConstraintRegion::outside:
+    case ConstraintRegion::inside:
+      return;
+  }
+  throw std::invalid_argument("checkpoint cylinder uses an unknown allowed region");
+}
+
 void validate_constraint_state(ConstraintId next_id, std::span<const PlaneConstraint> planes,
-                               std::span<const SphereConstraint> spheres) {
+                               std::span<const SphereConstraint> spheres,
+                               std::span<const BoxConstraint> boxes,
+                               std::span<const CylinderConstraint> cylinders) {
   if (next_id == invalid_constraint_id) {
     throw std::invalid_argument("checkpoint next constraint identifier is invalid");
   }
-  if (spheres.size() > std::numeric_limits<std::size_t>::max() - planes.size()) {
-    throw std::overflow_error("checkpoint constraint count overflow");
-  }
-  std::unordered_set<ConstraintId> ids;
-  ids.reserve(planes.size() + spheres.size());
-  ConstraintId previous_plane = invalid_constraint_id;
-  for (const auto& plane : planes) {
-    validate_plane(plane);
-    if (plane.id <= previous_plane || plane.id >= next_id) {
-      throw std::invalid_argument("checkpoint plane identifiers are not ordered and allocated");
-    }
-    if (!ids.insert(plane.id).second) {
-      throw std::invalid_argument("checkpoint contains a duplicate constraint identifier");
+  std::size_t total = planes.size();
+  for (const auto count : {spheres.size(), boxes.size(), cylinders.size()}) {
+    if (count > std::numeric_limits<std::size_t>::max() - total) {
+      throw std::overflow_error("checkpoint constraint count overflow");
     }
-    previous_plane = plane.id;
+    total += count;
   }
-  ConstraintId previous_sphere = invalid_constraint_id;
-  for (const auto& sphere : spheres) {
-    validate_sphere(sphere);
-    if (sphere.id <= previous_sphere || sphere.id >= next_id) {
-      throw std::invalid_argument("checkpoint sphere identifiers are not ordered and allocated");
-    }
-    if (!ids.insert(sphere.id).second) {
-      throw std::invalid_argument("checkpoint contains a duplicate constraint identifier");
+  std::unordered_set<ConstraintId> ids;
+  ids.reserve(total);
+  const auto check_ordered = [&ids, next_id](const auto& constraints, auto&& validate,
+                                             const char* kind) {
+    ConstraintId previous = invalid_constraint_id;
+    for (const auto& constraint : constraints) {
+      validate(constraint);
+      if (constraint.id <= previous || constraint.id >= next_id) {
+        throw std::invalid_argument(std::string("checkpoint ") + kind +
+                                    " identifiers are not ordered and allocated");
+      }
+      if (!ids.insert(constraint.id).second) {
+        throw std::invalid_argument("checkpoint contains a duplicate constraint identifier");
+      }
+      previous = constraint.id;
     }
-    previous_sphere = sphere.id;
-  }
+  };
+  check_ordered(planes, validate_plane, "plane");
+  check_ordered(spheres, validate_sphere, "sphere");
+  check_ordered(boxes, validate_box, "box");
+  check_ordered(cylinders, validate_cylinder, "cylinder");
 }
 
 std::size_t checked_offset_count(std::size_t cell_count) {
@@ -105,11 +143,15 @@ void validate_contact(const ExternalContact& contact, std::size_t cell_count) {
 }  // namespace
 
 void ConstraintSetCheckpoint::validate() const {
-  validate_constraint_state(next_id, planes, spheres);
+  validate_constraint_state(next_id, planes, spheres, boxes, cylinders);
 }
 
 ConstraintSet::ConstraintSet(const ConstraintSetCheckpoint& checkpoint)
-    : next_id_(checkpoint.next_id), planes_(checkpoint.planes), spheres_(checkpoint.spheres) {
+    : next_id_(checkpoint.next_id),
+      planes_(checkpoint.planes),
+      spheres_(checkpoint.spheres),
+      boxes_(checkpoint.boxes),
+      cylinders_(checkpoint.cylinders) {
   checkpoint.validate();
 }
 
@@ -155,25 +197,74 @@ ConstraintId ConstraintSet::add_sphere(const SphereConstraintInit& sphere) {
   return id;
 }
 
-std::size_t ConstraintSet::size() const noexcept { return planes_.size() + spheres_.size(); }
+ConstraintId ConstraintSet::add_box(const BoxConstraintInit& box) {
+  if (!finite(box.center) || !positive_finite_extents(box.half_extents)) {
+    throw std::invalid_argument("box geometry must be finite with positive half extents");
+  }
+  validate_coefficient(box.coefficient);
+  const auto id = allocate_id();
+  boxes_.push_back({
+      .id = id,
+      .center = box.center,
+      .half_extents = box.half_extents,
+      .coefficient = box.coefficient,
+      .allowed_region = box.allowed_region,
+  });
+  return id;
+}
 
-bool ConstraintSet::empty() const noexcept { return planes_.empty() && spheres_.empty(); }
+ConstraintId ConstraintSet::add_cylinder(const CylinderConstraintInit& cylinder) {
+  if (!finite(cylinder.center) || !std::isfinite(cylinder.radius) || cylinder.radius <= 0.0F ||
+      !std::isfinite(cylinder.half_height) || cylinder.half_height <= 0.0F) {
+    throw std::invalid_argument(
+        "cylinder geometry must be finite with a positive radius and half height");
+  }
+  validate_coefficient(cylinder.coefficient);
+  const auto id = allocate_id();
+  cylinders_.push_back({
+      .id = id,
+      .center = cylinder.center,
+      .radius = cylinder.radius,
+      .half_height = cylinder.half_height,
+      .coefficient = cylinder.coefficient,
+      .allowed_region = cylinder.allowed_region,
+  });
+  return id;
+}
+
+std::size_t ConstraintSet::size() const noexcept {
+  return planes_.size() + spheres_.size() + boxes_.size() + cylinders_.size();
+}
+
+bool ConstraintSet::empty() const noexcept {
+  return planes_.empty() && spheres_.empty() && boxes_.empty() && cylinders_.empty();
+}
 
 std::span<const PlaneConstraint> ConstraintSet::planes() const& noexcept { return planes_; }
 
 std::span<const SphereConstraint> ConstraintSet::spheres() const& noexcept { return spheres_; }
 
+std::span<const BoxConstraint> ConstraintSet::boxes() const& noexcept { return boxes_; }
+
+std::span<const CylinderConstraint> ConstraintSet::cylinders() const& noexcept {
+  return cylinders_;
+}
+
 ConstraintSetCheckpoint ConstraintSet::checkpoint() const {
   ConstraintSetCheckpoint result{
       .next_id = next_id_,
       .planes = planes_,
       .spheres = spheres_,
+      .boxes = boxes_,
+      .cylinders = cylinders_,
   };
   result.validate();
   return result;
 }
 
-void ConstraintSet::validate() const { validate_constraint_state(next_id_, planes_, spheres_); }
+void ConstraintSet::validate() const {
+  validate_constraint_state(next_id_, planes_, spheres_, boxes_, cylinders_);
+}
 
 void validate_constraint_contact_parameters(const ConstraintContactParameters& parameters) {
   if (!std::isfinite(parameters.activation_margin) || parameters.activation_margin < 0.0F) {

cpp/core/contact_graph.cpp

@@ -70,8 +70,8 @@ void validate_contact_parameters(const ContactParameters& parameters) {
   }
 }
 
-std::vector<ContactCandidate> find_cell_contact_candidates(
-    const WorldState& state, const ContactParameters& parameters) {
+std::vector<ContactCandidate> find_cell_contact_candidates(const WorldState& state,
+                                                           const ContactParameters& parameters) {
   validate_contact_parameters(parameters);
   const auto geometry = state.geometry_state();
   std::vector<CapsuleBounds> bounds;
@@ -112,10 +112,10 @@ std::vector<ContactCandidate> find_cell_contact_candidates(
     });
     active.erase(expired.begin(), expired.end());
     for (const auto* candidate : active) {
-      const auto overlaps_y = candidate->maximum_y >= current.minimum_y &&
-                              current.maximum_y >= candidate->minimum_y;
-      const auto overlaps_z = candidate->maximum_z >= current.minimum_z &&
-                              current.maximum_z >= candidate->minimum_z;
+      const auto overlaps_y =
+          candidate->maximum_y >= current.minimum_y && current.maximum_y >= candidate->minimum_y;
+      const auto overlaps_z =
+          candidate->maximum_z >= current.minimum_z && current.maximum_z >= candidate->minimum_z;
       if (!overlaps_y || !overlaps_z) {
         continue;
       }
@@ -125,11 +125,11 @@ std::vector<ContactCandidate> find_cell_contact_candidates(
     }
     active.push_back(&current);
   }
-  std::ranges::sort(candidates, [&geometry](const ContactCandidate& left,
-                                            const ContactCandidate& right) {
-    return std::tuple{geometry.ids[left.first_slot], geometry.ids[left.second_slot]} <
-           std::tuple{geometry.ids[right.first_slot], geometry.ids[right.second_slot]};
-  });
+  std::ranges::sort(
+      candidates, [&geometry](const ContactCandidate& left, const ContactCandidate& right) {
+        return std::tuple{geometry.ids[left.first_slot], geometry.ids[left.second_slot]} <
+               std::tuple{geometry.ids[right.first_slot], geometry.ids[right.second_slot]};
+      });
   return candidates;
 }
 

cpp/core/mechanics_integration.cpp

@@ -20,6 +20,8 @@ bool finite(const Vec3& value) {
   return std::isfinite(value.x) && std::isfinite(value.y) && std::isfinite(value.z);
 }
 
+}  // namespace
+
 Vec3 rotate_axis_angle(Vec3 direction, Vec3 rotation, float max_rotation) {
   const auto magnitude = norm(rotation);
   if (magnitude <= 1.0e-12F || max_rotation == 0.0F) {
@@ -33,8 +35,6 @@ Vec3 rotate_axis_angle(Vec3 direction, Vec3 rotation, float max_rotation) {
                     axis * (dot(axis, direction) * (1.0F - cosine)));
 }
 
-}  // namespace
-
 void validate_mechanics_integration_parameters(const MechanicsIntegrationParameters& parameters) {
   if (!std::isfinite(parameters.max_rotation_radians) || parameters.max_rotation_radians < 0.0F) {
     throw std::invalid_argument("mechanics rotation limit must be finite and non-negative");
@@ -74,10 +74,9 @@ void integrate_mechanics_result(WorldState& state, const MechanicsSolveResult& r
     const auto applied_length_increment =
         cell.fixed ? desired_increment : std::max(0.0F, desired_increment + correction.length);
     const auto new_position = cell.fixed ? cell.position : cell.position + correction.translation;
-    const auto new_direction =
-        cell.fixed ? cell.direction
-                   : rotate_axis_angle(cell.direction, correction.rotation,
-                                       parameters.max_rotation_radians);
+    const auto new_direction = cell.fixed ? cell.direction
+                                          : rotate_axis_angle(cell.direction, correction.rotation,
+                                                              parameters.max_rotation_radians);
     const auto new_length = cell.length + applied_length_increment;
     if (!finite(new_position) || !finite(new_direction) || !std::isfinite(new_length)) {
       throw std::overflow_error("mechanics integration produced non-finite geometry");