Skip to content

Commit ab52b57

Browse files
francoisbonneaufrancoisbonneau
committed
review arnaud
1 parent 9d53c3b commit ab52b57

5 files changed

Lines changed: 209 additions & 159 deletions

File tree

‎include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp‎

Lines changed: 80 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include <geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp>
3030
#include <geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp>
3131
#include <geode/stochastic/sampling/mcmc/models/components/density_term.hpp>
32+
#include <geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp>
3233
#include <geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp>
3334
#include <geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp>
3435
#include <geode/stochastic/spatial/pairwise_interactions.hpp>
@@ -107,7 +108,8 @@ namespace geode
107108
std::make_unique< ProposalKernel< OwnerSegment2D > >();
108109

109110
// Mapping set names -> UUID
110-
std::unordered_map< std::string, uuid > name_to_uuid;
111+
// std::unordered_map< std::string, uuid >
112+
// set_name_to_uuid_;
111113

112114
// Step 1: create object sets and samplers
113115
for( const auto& set_desc : set_descriptors_ )
@@ -120,7 +122,7 @@ namespace geode
120122
fixed_object[0], fixed_object[1] },
121123
set_id, true );
122124
}
123-
name_to_uuid[set_desc.name] = set_id;
125+
set_name_to_uuid_[set_desc.name] = set_id;
124126

125127
auto length_distribution =
126128
DoubleSampler::create_distribution( set_desc.length );
@@ -139,63 +141,12 @@ namespace geode
139141
// Step 2: create density energy terms
140142
for( const auto& set_desc : set_descriptors_ )
141143
{
142-
const auto set_id = name_to_uuid.at( set_desc.name );
143-
// p20
144-
this->ordered_energy_terms_.push_back(
145-
this->energy_terms_collection_.add_energy_term(
146-
std::make_unique< DensityTerm< OwnerSegment2D > >(
147-
absl::StrCat( set_desc.name, "_density" ),
148-
set_desc.p20, std::vector< uuid >{ set_id },
149-
this->domain_ ) ) );
150-
// intentisty
151-
if( std::fabs( set_desc.p21 - 1. ) > GLOBAL_EPSILON )
152-
{
153-
this->ordered_energy_terms_.push_back(
154-
this->energy_terms_collection_.add_energy_term(
155-
std::make_unique< IntensityTerm >(
156-
absl::StrCat( set_desc.name, "_intensity" ),
157-
set_desc.p21, std::vector< uuid >{ set_id },
158-
0.5 * this->domain_.smallest_length(),
159-
this->domain_ ) ) );
160-
}
161-
// spacing
162-
if( set_desc.minimal_spacing > GLOBAL_EPSILON )
163-
{
164-
auto interaction = std::make_unique<
165-
EuclideanCutoffInteraction< OwnerSegment2D > >(
166-
set_desc.minimal_spacing,
167-
PairwiseInteraction<
168-
OwnerSegment2D >::SCOPE::same_set );
169-
170-
this->ordered_energy_terms_.push_back(
171-
this->energy_terms_collection_.add_energy_term(
172-
std::make_unique< PairwiseTerm< OwnerSegment2D > >(
173-
absl::StrCat( set_desc.name, "_min_spacing" ),
174-
0., std::vector< uuid >{ set_id },
175-
std::move( interaction ), this->domain_ ) ) );
176-
}
144+
set_density_term( set_desc );
145+
set_intensity_term( set_desc );
146+
set_minimal_spacing_term( set_desc );
177147
}
178-
// x node monitoring
179-
if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON )
180-
{
181-
std::vector< uuid > set_uuids;
182-
set_uuids.reserve( name_to_uuid.size() );
183-
for( const auto& [name, id] : name_to_uuid )
184-
{
185-
set_uuids.push_back( id );
186-
}
187-
auto interaction = std::make_unique<
188-
EuclideanCutoffInteraction< OwnerSegment2D > >(
189-
0., PairwiseInteraction<
190-
OwnerSegment2D >::SCOPE::different_set );
148+
set_x_intersection_term();
191149

192-
this->ordered_energy_terms_.push_back(
193-
this->energy_terms_collection_.add_energy_term(
194-
std::make_unique< PairwiseTerm< OwnerSegment2D > >(
195-
absl::StrCat( "inter_set_x_nodes" ), beta_x_node_,
196-
set_uuids, std::move( interaction ),
197-
this->domain_ ) ) );
198-
}
199150
this->mh_sampler_ =
200151
std::make_unique< MetropolisHastings< OwnerSegment2D > >(
201152
this->energy_terms_collection_,
@@ -240,8 +191,80 @@ namespace geode
240191
return message;
241192
}
242193

194+
private:
195+
void set_density_term( const FractureSetDescription& set_desc )
196+
{
197+
const auto set_id = set_name_to_uuid_.at( set_desc.name );
198+
this->ordered_energy_terms_.push_back(
199+
this->energy_terms_collection_.add_energy_term(
200+
std::make_unique< DensityTerm< OwnerSegment2D > >(
201+
absl::StrCat( set_desc.name, "_density" ), set_desc.p20,
202+
std::vector< uuid >{ set_id }, this->domain_ ) ) );
203+
}
204+
205+
void set_intensity_term( const FractureSetDescription& set_desc )
206+
{
207+
if( std::fabs( set_desc.p21 - 1. ) < GLOBAL_EPSILON )
208+
{
209+
return;
210+
}
211+
const auto set_id = set_name_to_uuid_.at( set_desc.name );
212+
this->ordered_energy_terms_.push_back(
213+
this->energy_terms_collection_.add_energy_term(
214+
std::make_unique< IntensityTerm >(
215+
absl::StrCat( set_desc.name, "_intensity" ),
216+
set_desc.p21, std::vector< uuid >{ set_id },
217+
0.5 * this->domain_.smallest_length(),
218+
this->domain_ ) ) );
219+
}
220+
221+
void set_minimal_spacing_term( const FractureSetDescription& set_desc )
222+
{
223+
if( set_desc.minimal_spacing < GLOBAL_EPSILON )
224+
{
225+
return;
226+
}
227+
const auto set_id = set_name_to_uuid_.at( set_desc.name );
228+
auto interaction = std::make_unique<
229+
EuclideanCutoffInteraction< OwnerSegment2D > >(
230+
set_desc.minimal_spacing,
231+
PairwiseInteraction< OwnerSegment2D >::SCOPE::same_set );
232+
233+
this->ordered_energy_terms_.push_back(
234+
this->energy_terms_collection_.add_energy_term(
235+
std::make_unique< PairwiseTerm< OwnerSegment2D > >(
236+
absl::StrCat( set_desc.name, "_min_spacing" ), 0.,
237+
std::vector< uuid >{ set_id }, std::move( interaction ),
238+
this->domain_ ) ) );
239+
}
240+
241+
void set_x_intersection_term()
242+
{
243+
if( std::fabs( beta_x_node_ - 1. ) > GLOBAL_EPSILON )
244+
{
245+
std::vector< uuid > set_uuids;
246+
set_uuids.reserve( set_name_to_uuid_.size() );
247+
for( const auto& [name, id] : set_name_to_uuid_ )
248+
{
249+
set_uuids.push_back( id );
250+
}
251+
auto interaction = std::make_unique<
252+
EuclideanCutoffInteraction< OwnerSegment2D > >(
253+
0., PairwiseInteraction<
254+
OwnerSegment2D >::SCOPE::different_set );
255+
256+
this->ordered_energy_terms_.push_back(
257+
this->energy_terms_collection_.add_energy_term(
258+
std::make_unique< PairwiseTerm< OwnerSegment2D > >(
259+
absl::StrCat( "inter_set_x_nodes" ), beta_x_node_,
260+
set_uuids, std::move( interaction ),
261+
this->domain_ ) ) );
262+
}
263+
}
264+
243265
private:
244266
std::vector< FractureSetDescription > set_descriptors_;
267+
std::unordered_map< std::string, uuid > set_name_to_uuid_;
245268
// x node monitoring
246269
double beta_x_node_{ 1. };
247270
};

‎include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp‎

Lines changed: 0 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -26,59 +26,6 @@
2626

2727
#include <geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp>
2828

29-
double length_inside_box( const geode::Point2D& p0,
30-
const geode::Point2D& p1,
31-
const geode::BoundingBox2D& box )
32-
{
33-
double dx = p1.value( 0 ) - p0.value( 0 );
34-
double dy = p1.value( 1 ) - p0.value( 1 );
35-
36-
double t_min = 0.0;
37-
double t_max = 1.0;
38-
39-
auto update_interval = [&]( double p, double q_min, double q_max ) -> bool {
40-
if( p == 0.0 )
41-
{
42-
// Segment parallel to axis
43-
return ( q_min <= 0.0 && 0.0 <= q_max );
44-
}
45-
46-
double t1 = q_min / p;
47-
double t2 = q_max / p;
48-
if( t1 > t2 )
49-
std::swap( t1, t2 );
50-
51-
t_min = std::max( t_min, t1 );
52-
t_max = std::min( t_max, t2 );
53-
54-
return t_min <= t_max;
55-
};
56-
57-
// X axis
58-
if( !update_interval( dx, box.min().value( 0 ) - p0.value( 0 ),
59-
box.max().value( 0 ) - p0.value( 0 ) ) )
60-
{
61-
return 0.0;
62-
}
63-
64-
// Y axis
65-
if( !update_interval( dy, box.min().value( 1 ) - p0.value( 1 ),
66-
box.max().value( 1 ) - p0.value( 1 ) ) )
67-
{
68-
return 0.0;
69-
}
70-
71-
if( t_max <= t_min )
72-
{
73-
return 0.0;
74-
}
75-
76-
double clipped_dx = dx * ( t_max - t_min );
77-
double clipped_dy = dy * ( t_max - t_min );
78-
79-
return std::sqrt( clipped_dx * clipped_dx + clipped_dy * clipped_dy );
80-
}
81-
8229
namespace geode
8330
{
8431
template < typename ObjectType >
@@ -112,33 +59,4 @@ namespace geode
11259
{
11360
}
11461
};
115-
116-
class IntensityTerm
117-
: public SingleObjectTerm< OwnerSegment2D,
118-
std::function< double( const OwnerSegment2D&,
119-
const SpatialDomain< OwnerSegment2D::dim >& ) > >
120-
{
121-
public:
122-
explicit IntensityTerm( std::string_view name,
123-
double lambda,
124-
std::vector< uuid > targeted_set_ids,
125-
double caracteristic_length,
126-
const SpatialDomain< OwnerSegment2D::dim >& domain )
127-
: SingleObjectTerm< OwnerSegment2D,
128-
std::function< double( const OwnerSegment2D&,
129-
const SpatialDomain< OwnerSegment2D::dim >& ) > >(
130-
name,
131-
lambda,
132-
std::move( targeted_set_ids ),
133-
1.0 / caracteristic_length,
134-
[]( const OwnerSegment2D& segment,
135-
const SpatialDomain< OwnerSegment2D::dim >& domain ) {
136-
auto seg_extremities = segment.vertices();
137-
return length_inside_box( seg_extremities[0],
138-
seg_extremities[1], domain.box() );
139-
},
140-
domain )
141-
{
142-
}
143-
};
14462
} // namespace geode

‎include/geode/stochastic/sampling/mcmc/models/components/intensity_term.hpp‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
/*
2+
* Copyright (c) 2019 - 2025 Geode-solutions
3+
*
4+
* Permission is hereby granted, free of charge, to any person obtaining a copy
5+
* of this software and associated documentation files (the "Software"), to deal
6+
* in the Software without restriction, including without limitation the rights
7+
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8+
* copies of the Software, and to permit persons to whom the Software is
9+
* furnished to do so, subject to the following conditions:
10+
*
11+
* The above copyright notice and this permission notice shall be included in
12+
* all copies or substantial portions of the Software.
13+
*
14+
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15+
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16+
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17+
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18+
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19+
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20+
* SOFTWARE.
21+
*
22+
*/
23+
#pragma once
24+
25+
#include <geode/stochastic/spatial/object_sets.hpp>
26+
27+
#include <geode/stochastic/sampling/mcmc/models/components/single_object_term.hpp>
28+
29+
double length_inside_box( const geode::Point2D& segment_start,
30+
const geode::Point2D& segment_end,
31+
const geode::BoundingBox2D& box )
32+
{
33+
// Segment direction
34+
const double dir_x = segment_end.value( 0 ) - segment_start.value( 0 );
35+
const double dir_y = segment_end.value( 1 ) - segment_start.value( 1 );
36+
37+
// Parameter interval where the segment is inside the box
38+
double t_enter = 0.0; // start of valid interval
39+
double t_exit = 1.0; // end of valid interval
40+
41+
// Clips the parametric segment against one axis interval
42+
auto clip_against_axis_interval = [&t_enter, &t_exit]( double direction,
43+
double min_distance,
44+
double max_distance ) -> bool {
45+
// Segment is parallel to this axis
46+
if( std::fabs( direction ) < geode::GLOBAL_EPSILON )
47+
{
48+
// Outside the interval → no intersection
49+
return min_distance <= 0.0 && 0.0 <= max_distance;
50+
}
51+
52+
double t0 = min_distance / direction;
53+
double t1 = max_distance / direction;
54+
55+
if( t0 > t1 )
56+
{
57+
std::swap( t0, t1 );
58+
}
59+
60+
t_enter = std::max( t_enter, t0 );
61+
t_exit = std::min( t_exit, t1 );
62+
63+
return t_enter <= t_exit;
64+
};
65+
66+
// Clip against X interval of the box
67+
if( !clip_against_axis_interval( dir_x,
68+
box.min().value( 0 ) - segment_start.value( 0 ),
69+
box.max().value( 0 ) - segment_start.value( 0 ) ) )
70+
{
71+
return 0.0;
72+
}
73+
74+
// Clip against Y interval of the box
75+
if( !clip_against_axis_interval( dir_y,
76+
box.min().value( 1 ) - segment_start.value( 1 ),
77+
box.max().value( 1 ) - segment_start.value( 1 ) ) )
78+
{
79+
return 0.0;
80+
}
81+
82+
// No portion of the segment is inside the box
83+
if( t_exit <= t_enter )
84+
{
85+
return 0.0;
86+
}
87+
88+
// Length of the clipped segment
89+
const double clipped_dx = dir_x * ( t_exit - t_enter );
90+
const double clipped_dy = dir_y * ( t_exit - t_enter );
91+
92+
return std::sqrt( clipped_dx * clipped_dx + clipped_dy * clipped_dy );
93+
}
94+
95+
namespace geode
96+
{
97+
class IntensityTerm
98+
: public SingleObjectTerm< OwnerSegment2D,
99+
std::function< double( const OwnerSegment2D&,
100+
const SpatialDomain< OwnerSegment2D::dim >& ) > >
101+
{
102+
public:
103+
explicit IntensityTerm( std::string_view name,
104+
double lambda,
105+
std::vector< uuid > targeted_set_ids,
106+
double caracteristic_length,
107+
const SpatialDomain< OwnerSegment2D::dim >& domain )
108+
: SingleObjectTerm< OwnerSegment2D,
109+
std::function< double( const OwnerSegment2D&,
110+
const SpatialDomain< OwnerSegment2D::dim >& ) > >(
111+
name,
112+
lambda,
113+
std::move( targeted_set_ids ),
114+
1.0 / caracteristic_length,
115+
[]( const OwnerSegment2D& segment,
116+
const SpatialDomain< OwnerSegment2D::dim >& domain ) {
117+
auto seg_extremities = segment.vertices();
118+
return length_inside_box( seg_extremities[0],
119+
seg_extremities[1], domain.box() );
120+
},
121+
domain )
122+
{
123+
}
124+
};
125+
} // namespace geode

0 commit comments

Comments
 (0)