feat(Helpers): Fracture Simulation Helpers

include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp

@@ -81,7 +81,7 @@ namespace geode
         }
 
     private:
-        const BoundingBox< dimension >& box_;
+        BoundingBox< dimension > box_;
     };
 
 } // namespace geode

include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp

@@ -0,0 +1,103 @@
+/*
+ * Copyright (c) 2019 - 2025 Geode-solutions
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ */
+
+#pragma once
+
+#include <geode/geometry/basic_objects/segment.hpp>
+#include <geode/geometry/basic_objects/sphere.hpp>
+
+#include <geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp>
+#include <geode/stochastic/sampling/direct/point_uniform_sampler.hpp>
+#include <geode/stochastic/sampling/direct/segment_uniform_sampler.hpp>
+
+namespace geode
+{
+    class UniformSegmentSetSampler : public ObjectSetSampler< OwnerSegment2D >
+    {
+    public:
+        UniformSegmentSetSampler( const BoundingBox< 2 >& box,
+            const DoubleSampler::Distribution& length,
+            const DoubleSampler::Distribution& azimuth )
+            : ObjectSetSampler< OwnerSegment2D >{},
+              box_{ box },
+              length_{ length },
+              azimuth_{ azimuth }
+        {
+            auto volume = box_.n_volume();
+            OPENGEODE_EXCEPTION( volume != 0.,
+                "[SegmentSetSampler] - Undefined Bounding Box (volume == 0)." );
+            this->log_pdf_ = -std::log( volume );
+        }
+
+        OwnerSegment2D sample( RandomEngine& engine ) const override
+        {
+            auto seg = SegmentUniformSampler::sample(
+                engine, box_, length_, azimuth_ );
+            return seg;
+        }
+
+        OwnerSegment2D change(
+            const OwnerSegment2D& obj, RandomEngine& engine ) const override
+        {
+            double ratio = 0.1;
+            const auto& extremities = obj.vertices();
+            const auto current =
+                static_cast< local_index_t >( engine.sample_bernoulli( 0.5 ) );
+
+            geode::Sphere< 2 > ball{ extremities[current],
+                ratio * obj.length() };
+
+            auto new_point = PointUniformSampler::sample< 2 >( engine, ball );
+            constexpr index_t max_try{ 100 };
+            for( const auto try_id : geode::Range{ max_try } )
+            {
+                if( box_.contains( new_point ) )
+                {
+                    OwnerSegment2D new_segment{ obj };
+                    new_segment.set_point( current, new_point );
+                    return new_segment;
+                }
+                new_point = PointUniformSampler::sample< 2 >( engine, ball );
+            }
+            throw OpenGeodeException( absl::StrCat(
+                "[PointSampler] - Cannot find a point in the box: ",
+                box_.string() ) );
+            return obj;
+        }
+
+    private:
+        bool is_valid_object( const OwnerSegment2D& obj ) const override
+        {
+            const auto& extremities = obj.vertices();
+
+            return box_.contains( extremities[0] )
+                   && box_.contains( extremities[1] );
+        }
+
+    private:
+        BoundingBox2D box_;
+        DoubleSampler::Distribution length_;
+        DoubleSampler::Distribution azimuth_;
+    };
+
+} // namespace geode

include/geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp

@@ -219,8 +219,6 @@ namespace geode
             const auto old_object_id = proposal.old_object_id();
             const auto delta_log_energy = gibbs_energy_.delta_log_change(
                 state, old_object_id, new_object );
-            // should we test that objects are in the same group?
-            // should be ensured by the dynamic
             return accept_or_reject( proposal, state, engine, delta_log_energy,
                 []( auto& state, auto& proposal ) {
                     state.update_object( proposal.old_object_id(),

include/geode/stochastic/spatial/pairwise_interactions.hpp

@@ -25,6 +25,7 @@
 
 #include <geode/basic/uuid.hpp>
 
+#include <geode/geometry/basic_objects/segment.hpp>
 #include <geode/geometry/distance.hpp>
 #include <geode/geometry/point.hpp>
 
@@ -74,38 +75,26 @@ namespace geode
         SCOPE scope_{ SCOPE::all_set };
     };
 
+    /*!
+     * EuclideanCutoffInteraction
+     * A pairwise interaction that returns 1 if the Euclidean distance
+     * between objects is within a cutoff radius, otherwise 0.
+     */
     template < typename Type >
     class EuclideanCutoffInteraction : public PairwiseInteraction< Type >
     {
     public:
-        EuclideanCutoffInteraction( double cutoff_distance )
-            : PairwiseInteraction< Type >(), cutoff_distance_( cutoff_distance )
-        {
-        }
+        explicit EuclideanCutoffInteraction( double cutoff_distance );
         EuclideanCutoffInteraction( double cutoff_distance,
-            typename PairwiseInteraction< Type >::SCOPE scope )
-            : PairwiseInteraction< Type >( scope ),
-              cutoff_distance_( cutoff_distance )
-        {
-        }
+            typename PairwiseInteraction< Type >::SCOPE scope );
 
-        double neighborhood_searching_distance() const override
-        {
-            return cutoff_distance_;
-        }
+        double neighborhood_searching_distance() const override;
 
     protected:
         double compute( const ObjectRef< Type >& object_a,
-            const ObjectRef< Type >& object_b ) const override
-        {
-            auto dist = geode::point_point_distance(
-                object_barycenter( object_a.object ),
-                object_barycenter( object_b.object ) );
-            return dist <= cutoff_distance_ ? 1.0 : 0.0;
-        }
+            const ObjectRef< Type >& object_b ) const override;
 
     private:
         double cutoff_distance_{ GLOBAL_EPSILON };
     };
-
 } // namespace geode

src/geode/stochastic/CMakeLists.txt

@@ -25,6 +25,7 @@ add_geode_library(
         "spatial/object_set.cpp"
         "spatial/object_sets.cpp"
         "spatial/object_neighborhood.cpp"
+        "spatial/pairwise_interactions.cpp"
         "sampling/direct/ball_sampler.cpp"
         "sampling/direct/bounding_box_sampler.cpp"
         "sampling/direct/double_sampler.cpp"
@@ -43,6 +44,7 @@ add_geode_library(
         "spatial/details/RTree.hpp"
         "sampling/direct/object_set_sampler/object_set_sampler.hpp"
         "sampling/direct/object_set_sampler/point_set_sampler.hpp"
+        "sampling/direct/object_set_sampler/segment_set_sampler.hpp"
         "sampling/direct/ball_sampler.hpp"
         "sampling/direct/bounding_box_sampler.hpp"
         "sampling/direct/double_sampler.hpp"

src/geode/stochastic/spatial/object_set.cpp

@@ -87,8 +87,8 @@ namespace geode
     template < typename Type >
     std::string ObjectSet< Type >::string() const
     {
-        return absl::StrCat(
-            "ObjectSet ", this->id().string(), " with ", size(), " objects" );
+        return absl::StrCat( "ObjectSet ", this->name(), " (",
+            this->id().string(), ") contains ", size(), " objects" );
     }
 
     // Explicit template instantiation

src/geode/stochastic/spatial/object_sets.cpp

@@ -4,31 +4,7 @@
 
 #include <geode/geometry/basic_objects/segment.hpp>
 #include <geode/geometry/point.hpp>
-namespace
-{
-    template < geode::index_t dimension >
-    geode::Vector< dimension > enlarge_vector( double distance );
-
-    template <>
-    geode::Vector2D enlarge_vector< 2 >( double distance )
-    {
-        return geode::Vector2D{ { distance, distance } };
-    }
 
-    template <>
-    geode::Vector3D enlarge_vector< 3 >( double distance )
-    {
-        return geode::Vector3D{ { distance, distance, distance } };
-    }
-    template < geode::index_t dimension >
-    void enlarge_box(
-        geode::BoundingBox< dimension >& box, double search_distance )
-    {
-        auto vec = enlarge_vector< dimension >( search_distance );
-        box.add_point( box.min() - vec );
-        box.add_point( box.max() + vec );
-    }
-} // namespace
 namespace geode
 {
     template < typename Type >
@@ -150,7 +126,7 @@ namespace geode
         const ObjectId& object_id, double searching_distance ) const
     {
         auto box = object_bounding_box( get_object( object_id ) );
-        enlarge_box< Type::dim >( box, searching_distance );
+        box.extends( searching_distance * 2. );
         return neighborhood_.get_all_neighbor_ids( box, object_id );
     }
 
@@ -159,7 +135,7 @@ namespace geode
         const Type& object, double searching_distance ) const
     {
         auto box = object_bounding_box( object );
-        enlarge_box< Type::dim >( box, searching_distance );
+        box.extends( searching_distance * 2. );
         return neighborhood_.get_all_neighbor_ids( box, std::nullopt );
     }
 
@@ -173,12 +149,11 @@ namespace geode
     template < typename Type >
     std::string ObjectSets< Type >::string() const
     {
-        auto message =
-            absl::StrCat( "ObjectSets with ", nb_sets(), " ObjectSet:" );
+        auto message = absl::StrCat( "ObjectSets with ", nb_objects(),
+            " objects in, ", nb_sets(), " sets" );
         for( const auto& [set_id, objs] : sets_ )
         {
-            absl::StrAppend( &message, "\n\t - set uuid: ", set_id.string(),
-                "--> ", objs.string() );
+            absl::StrAppend( &message, "\n\t --> ", objs.string() );
         }
         return message;
     }

src/geode/stochastic/spatial/pairwise_interactions.cpp

@@ -0,0 +1,86 @@
+/*
+ * Copyright (c) 2019 - 2025 Geode-solutions
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ */
+
+#include <geode/stochastic/spatial/pairwise_interactions.hpp>
+
+#include <geode/geometry/basic_objects/segment.hpp>
+#include <geode/geometry/distance.hpp>
+
+namespace
+{
+    template < geode::index_t dimension >
+    double compute_euclidean_distance( const geode::Point< dimension >& point1,
+        const geode::Point< dimension >& point2 )
+    {
+        return geode::point_point_distance( point1, point2 );
+    }
+    template < geode::index_t dimension >
+    double compute_euclidean_distance( const geode::Segment< dimension >& seg1,
+        const geode::Segment< dimension >& seg2 )
+    {
+        return std::get< 0 >( geode::segment_segment_distance( seg1, seg2 ) );
+    }
+} // namespace
+namespace geode
+{
+    template < typename Type >
+    EuclideanCutoffInteraction< Type >::EuclideanCutoffInteraction(
+        double cutoff_distance )
+        : PairwiseInteraction< Type >(), cutoff_distance_( cutoff_distance )
+    {
+    }
+
+    template < typename Type >
+    EuclideanCutoffInteraction< Type >::EuclideanCutoffInteraction(
+        double cutoff_distance,
+        typename PairwiseInteraction< Type >::SCOPE scope )
+        : PairwiseInteraction< Type >( scope ),
+          cutoff_distance_( cutoff_distance )
+    {
+    }
+
+    template < typename Type >
+    double EuclideanCutoffInteraction< Type >::neighborhood_searching_distance()
+        const
+    {
+        return cutoff_distance_;
+    }
+
+    template < typename Type >
+    double EuclideanCutoffInteraction< Type >::compute(
+        const ObjectRef< Type >& object_a,
+        const ObjectRef< Type >& object_b ) const
+    {
+        auto dist = compute_euclidean_distance< Type::dim >(
+            object_a.object, object_b.object );
+        return dist <= cutoff_distance_ ? 1.0 : 0.0;
+    }
+
+    template class opengeode_stochastic_stochastic_api
+        EuclideanCutoffInteraction< Point< 2 > >;
+    template class opengeode_stochastic_stochastic_api
+        EuclideanCutoffInteraction< Point< 3 > >;
+
+    template class opengeode_stochastic_stochastic_api
+        EuclideanCutoffInteraction< OwnerSegment< 2 > >;
+} // namespace geode

tests/stochastic/CMakeLists.txt

@@ -121,6 +121,14 @@ add_geode_test(
         ${PROJECT_NAME}::stochastic
 )
 
+add_geode_test(
+     SOURCE "sampling/mcmc/test-mh-fractures.cpp"
+     DEPENDENCIES
+         OpenGeode::basic
+         OpenGeode::geometry
+         ${PROJECT_NAME}::stochastic
+)
+
 add_geode_test(
      SOURCE "sampling/mcmc/test-mh-poisson.cpp"
      DEPENDENCIES

tests/stochastic/sampling/direct/test-double-sampler.cpp

@@ -51,7 +51,7 @@ void test_double_sampler( geode::RandomEngine& engine )
 
     geode::TruncatedGaussian t_gaussian_double;
     t_gaussian_double.mean = 0.;
-    t_gaussian_double.mean = 1.;
+    t_gaussian_double.standard_deviation = 1.;
     t_gaussian_double.max_value = 1.;
     t_gaussian_double.min_value = -1.;
     value = 1000.;