From 921cfdd9eb3295ecb8410dc257ff3651308485d2 Mon Sep 17 00:00:00 2001 From: nicogno Date: Fri, 2 Aug 2024 10:23:26 -0400 Subject: [PATCH 1/5] Switched absolute to relative additional radius --- src/core/interaction_force.cc | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/core/interaction_force.cc b/src/core/interaction_force.cc index b8c463d89..34c039199 100644 --- a/src/core/interaction_force.cc +++ b/src/core/interaction_force.cc @@ -62,21 +62,19 @@ void InteractionForce::ForceBetweenSpheres(const Agent* sphere_lhs, Real3* result) const { const Real3& ref_mass_location = sphere_lhs->GetPosition(); real_t ref_diameter = sphere_lhs->GetDiameter(); - real_t ref_iof_coefficient = 0.15; const Real3& nb_mass_location = sphere_rhs->GetPosition(); real_t nb_diameter = sphere_rhs->GetDiameter(); - real_t nb_iof_coefficient = 0.15; auto c1 = ref_mass_location; real_t r1 = 0.5 * ref_diameter; auto c2 = nb_mass_location; real_t r2 = 0.5 * nb_diameter; // We take virtual bigger radii to have a distant interaction, to get a - // desired density. - real_t additional_radius = - 10.0 * std::min(ref_iof_coefficient, nb_iof_coefficient); - r1 += additional_radius; - r2 += additional_radius; + // desired density. We assume an interaction distance 3% larger then the + // agent's radius + real_t additional_radius_scaling = 1.03; + r1 *= additional_radius_scaling; + r2 *= additional_radius_scaling; // the 3 components of the vector c2 -> c1 real_t comp1 = c1[0] - c2[0]; real_t comp2 = c1[1] - c2[1]; From 42f50176b5339ce82c9c95cfcd2762735200fdd9 Mon Sep 17 00:00:00 2001 From: nicogno Date: Tue, 6 Aug 2024 17:31:32 -0400 Subject: [PATCH 2/5] Updated test values --- test/unit/core/interaction_force_test.cc | 14 ++--- .../operation/mechanical_forces_op_test.h | 59 ++++++++++--------- 2 files changed, 37 insertions(+), 36 deletions(-) diff --git a/test/unit/core/interaction_force_test.cc b/test/unit/core/interaction_force_test.cc index 5ae3e45f1..de6efe2bc 100644 --- a/test/unit/core/interaction_force_test.cc +++ b/test/unit/core/interaction_force_test.cc @@ -37,20 +37,20 @@ TEST(InteractionForce, GeneralSphere) { InteractionForce force; auto result = force.Calculate(&cell, &nb); - EXPECT_NEAR(7.1429184067241138, result[0], abs_error::value); - EXPECT_NEAR(6.4935621879310119, result[1], abs_error::value); - EXPECT_NEAR(5.8442059691379109, result[2], abs_error::value); + EXPECT_NEAR(4.50152537445543206341835066295, result[0], abs_error::value); + EXPECT_NEAR(4.09229579495948369401668242086, result[1], abs_error::value); + EXPECT_NEAR(3.68306621546353532461501417877, result[2], abs_error::value); nb.SetDiameter(10); nb.SetPosition({5, 5, 0}); result = force.Calculate(&cell, &nb); - EXPECT_NEAR(-5.7454658831720176, result[0], abs_error::value); - EXPECT_NEAR(-5.892785521202069, result[1], abs_error::value); - EXPECT_NEAR(1.3258767422704656, result[2], abs_error::value); + EXPECT_NEAR(-2.99648347356091527350007596923, result[0], abs_error::value); + EXPECT_NEAR(-3.07331638313940028051289842998, result[1], abs_error::value); + EXPECT_NEAR(0.691496186206365063115402146746, result[2], abs_error::value); } -/// Tests the special case that non of the neighbors overlap +/// Tests the special case that none of the neighbors overlap /// with the reference cell TEST(InteractionForce, AllNonOverlappingSphere) { Simulation simulation(TEST_NAME); diff --git a/test/unit/core/operation/mechanical_forces_op_test.h b/test/unit/core/operation/mechanical_forces_op_test.h index 55b4efa14..a124f8ee7 100644 --- a/test/unit/core/operation/mechanical_forces_op_test.h +++ b/test/unit/core/operation/mechanical_forces_op_test.h @@ -59,13 +59,13 @@ inline void RunTest(const std::string& environment) { Cell* final_cell1 = dynamic_cast(rm->GetAgent(ref_uid + 1)); auto final_position = final_cell0->GetPosition(); EXPECT_NEAR(0, final_position[0], abs_error::value); - EXPECT_NEAR(-0.07797206232558615, final_position[1], + EXPECT_NEAR(-0.0494590384253325649715244643239, final_position[1], abs_error::value); EXPECT_NEAR(0, final_position[2], abs_error::value); // cell 1 final_position = final_cell1->GetPosition(); EXPECT_NEAR(0, final_position[0], abs_error::value); - EXPECT_NEAR(5.0980452768658333, final_position[1], abs_error::value); + EXPECT_NEAR(5.06220489350739084993450043382, final_position[1], abs_error::value); EXPECT_NEAR(0, final_position[2], abs_error::value); // check if tractor_force has been reset to zero @@ -129,33 +129,34 @@ inline void RunTest2(const std::string& environment) { } // clang-format off - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 0)->GetPosition(), {-0.20160966809506442, -0.20160966809506442, -0.20160966809506442}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 1)->GetPosition(), {19.996726567569418, -0.22266672163763598, -0.22266672163763598}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 2)->GetPosition(), {40.201498173927305, -0.19986951265938641, -0.19986951265938641}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 3)->GetPosition(), {-0.22107235667371566, 19.99536937981701, -0.22243298066161277}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 4)->GetPosition(), {19.996445694277536, 19.991084649064845, -0.24302072168423453}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 5)->GetPosition(), {40.220941141351886, 19.995420431150347, -0.22054684872504365}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 6)->GetPosition(), {-0.19982869510573351, 40.201437349791846, -0.19959785780667869}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 7)->GetPosition(), {19.99678197072755, 40.219229813316794, -0.22031963374181529}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 8)->GetPosition(), {40.199677082903456, 40.19967381378099, -0.19789135788064738}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 9)->GetPosition(), {-0.2208401632148318, -0.22084860029084805, 19.994023485540392}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 10)->GetPosition(), {19.996453518945092, -0.2413353410797266, 19.989660774799667}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 11)->GetPosition(), {40.220709535497228, -0.21898884436167229, 19.994096990730121}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 12)->GetPosition(), {-0.23959752549368377, 19.995004901057818, 19.989639428390195}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 13)->GetPosition(), {19.996198530091725, 19.99036483907971, 19.984534769039289}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 14)->GetPosition(), {40.239451672678122, 19.995058936152933, 19.989784932711256}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 15)->GetPosition(), {-0.21895170976356038, 40.220651199374309, 19.994102585380165}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 16)->GetPosition(), {19.996511705095408, 40.237587787915302, 19.989763131796483}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 17)->GetPosition(), {40.218778612717593, 40.218767244773701, 19.994175099676564}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 18)->GetPosition(), {-0.19955064743630646, -0.19955761561629054, 40.201369950961137}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 19)->GetPosition(), {19.99679050342159, -0.22028367546827032, 40.219170471740057}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 20)->GetPosition(), {40.199399442577274, -0.19785191118676104, 40.199607855810939}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 21)->GetPosition(), {-0.2186842147430082, 19.995465190535288, 40.218898427254011}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 22)->GetPosition(), {19.996519321027819, 19.991238784452595, 40.235673745336591}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 23)->GetPosition(), {40.218511916378013, 19.995514817696375, 40.217044778537549}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 24)->GetPosition(), {-0.19780590198157336, 40.19933398679354, 40.199323871317965}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 25)->GetPosition(), {19.996843796541615, 40.216799284843852, 40.216779234925141}); - EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 26)->GetPosition(), {40.197616806612238, 40.197607143403182, 40.197597121203316}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 0)->GetPosition(), {-0.136650834179647427449523566794, -0.136650834179647427449523566794, -0.136650834179647427449523566794}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 1)->GetPosition(), {19.9978236653235009314379819452, -0.141015306673735183678484988771, -0.141015306673735183678484988771}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 2)->GetPosition(), {40.1365944997387418724705928750, -0.13576125310507888313992208380, -0.13576125310507888313992208380}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 3)->GetPosition(), {-0.14020905712356290333724784449, 19.9970489222767306358984663823, -0.140984279029959657535431602871}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 4)->GetPosition(), {19.9977851373048793139637014758, 19.9947099058674891346242305069, -0.145277624139375911514663940516}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 5)->GetPosition(), {40.1401388185864714084643836652, 19.9970715249564464264518999951, -0.14007326379892509482257787044}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 6)->GetPosition(), {-0.13574486720051072606503005095, 40.1365627529316842850017320254, -0.135723902739067722287343804703}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 7)->GetPosition(), {19.9978515057522633535842662641, 40.1392617421792200453064082392, -0.140043337731872570780977206275}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 8)->GetPosition(), {40.1356651625421862223149400611, 40.1356589129132839546632435903, -0.134846375207706975815553682423}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 9)->GetPosition(), {-0.140178174501168932165811283971, -0.140182476858812768317650143102, 19.9962788234873976200766333855}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 10)->GetPosition(), {19.9977862308751720034084509786, -0.144471737809141948159807531113, 19.9939369302356357458986222506}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 11)->GetPosition(), {40.1401080571319587799692461333, -0.139281644960194857517683474124, 19.9963112337827742562800909788}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 12)->GetPosition(), {-0.143636884973078624410733274983, 19.9970003497588359473907353340, 19.9939216324133267598249705192}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 13)->GetPosition(), {19.9977509281001790744848931789, 19.9946123129189552912189860276, 19.9914755278624715229091673136}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 14)->GetPosition(), {40.1435533751884737842596265535, 19.9970230447662831556627447542, 19.9939839395317135001120125930}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 15)->GetPosition(), {-0.139264372505880738304090978511, 40.1400752445670986163454723481, 19.9963056313964950800704680887}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 16)->GetPosition(), {19.9978140032122536103189974908, 40.1426469666206381715941682442, 19.9939685327911046231094964962}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 17)->GetPosition(), {40.1391708994512008085005089423, 40.1391603959532489256165370155, 19.9963378247512290443404377838}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 18)->GetPosition(), {-0.135703460393067269278536235005, -0.135707572042464525088701330238, 40.1365270491528043523269029624}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 19)->GetPosition(), {19.9978526859614648754687984712, -0.140029334896406153010924586016, 40.1392258289244541720704728805}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 20)->GetPosition(), {40.1356238143556917587581200186, -0.13483030541778221502463362419, 40.1356238215903481357585777958}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 21)->GetPosition(), {-0.139216463603416178525167668456, 19.9970888470926410320886260158, 40.1391768764013402113023482055}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 22)->GetPosition(), {19.9978151330991056052596669472, 19.9947714402862555915866948948, 40.1417125181374631632487142307}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 23)->GetPosition(), {40.1391231591415114231657215943, 19.9971109750473197461909678817, 40.1382739675789673370998391272}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 24)->GetPosition(), {-0.134810200828153573232552824888, 40.1355886055729997034154133547, 40.1355781327590436272539263897}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 25)->GetPosition(), {19.9978797095407642059976814643, 40.1382410601235742697478311222, 40.1382260635760805291432435869}); + EXPECT_ARR_NEAR(rm->GetAgent(ref_uid + 26)->GetPosition(), {40.1347081112982584823528491696, 40.1346979778000455604454679675, 40.1346876663558284147859833527}); + // clang-format on delete mechanical_forces_op; From 7e4b89e2c4278bde2d384c4f5c14ab4218a228ca Mon Sep 17 00:00:00 2001 From: nicogno Date: Tue, 6 Aug 2024 17:35:32 -0400 Subject: [PATCH 3/5] Format --- test/unit/core/interaction_force_test.cc | 18 ++++++++++++------ .../core/operation/mechanical_forces_op_test.h | 3 ++- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/test/unit/core/interaction_force_test.cc b/test/unit/core/interaction_force_test.cc index de6efe2bc..71b9e2eec 100644 --- a/test/unit/core/interaction_force_test.cc +++ b/test/unit/core/interaction_force_test.cc @@ -37,17 +37,23 @@ TEST(InteractionForce, GeneralSphere) { InteractionForce force; auto result = force.Calculate(&cell, &nb); - EXPECT_NEAR(4.50152537445543206341835066295, result[0], abs_error::value); - EXPECT_NEAR(4.09229579495948369401668242086, result[1], abs_error::value); - EXPECT_NEAR(3.68306621546353532461501417877, result[2], abs_error::value); + EXPECT_NEAR(4.50152537445543206341835066295, result[0], + abs_error::value); + EXPECT_NEAR(4.09229579495948369401668242086, result[1], + abs_error::value); + EXPECT_NEAR(3.68306621546353532461501417877, result[2], + abs_error::value); nb.SetDiameter(10); nb.SetPosition({5, 5, 0}); result = force.Calculate(&cell, &nb); - EXPECT_NEAR(-2.99648347356091527350007596923, result[0], abs_error::value); - EXPECT_NEAR(-3.07331638313940028051289842998, result[1], abs_error::value); - EXPECT_NEAR(0.691496186206365063115402146746, result[2], abs_error::value); + EXPECT_NEAR(-2.99648347356091527350007596923, result[0], + abs_error::value); + EXPECT_NEAR(-3.07331638313940028051289842998, result[1], + abs_error::value); + EXPECT_NEAR(0.691496186206365063115402146746, result[2], + abs_error::value); } /// Tests the special case that none of the neighbors overlap diff --git a/test/unit/core/operation/mechanical_forces_op_test.h b/test/unit/core/operation/mechanical_forces_op_test.h index a124f8ee7..5d6e129f2 100644 --- a/test/unit/core/operation/mechanical_forces_op_test.h +++ b/test/unit/core/operation/mechanical_forces_op_test.h @@ -65,7 +65,8 @@ inline void RunTest(const std::string& environment) { // cell 1 final_position = final_cell1->GetPosition(); EXPECT_NEAR(0, final_position[0], abs_error::value); - EXPECT_NEAR(5.06220489350739084993450043382, final_position[1], abs_error::value); + EXPECT_NEAR(5.06220489350739084993450043382, final_position[1], + abs_error::value); EXPECT_NEAR(0, final_position[2], abs_error::value); // check if tractor_force has been reset to zero From 782b5d45045d5671e292c02a0c338086e5cdd9ab Mon Sep 17 00:00:00 2001 From: nicogno Date: Wed, 7 Aug 2024 19:10:51 -0400 Subject: [PATCH 4/5] Made virtual agent radius scale factor adjustable by user --- src/core/interaction_force.cc | 10 ++++++---- src/core/param/param.cc | 2 ++ src/core/param/param.h | 8 ++++++++ 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/core/interaction_force.cc b/src/core/interaction_force.cc index 34c039199..27418e44d 100644 --- a/src/core/interaction_force.cc +++ b/src/core/interaction_force.cc @@ -64,6 +64,8 @@ void InteractionForce::ForceBetweenSpheres(const Agent* sphere_lhs, real_t ref_diameter = sphere_lhs->GetDiameter(); const Real3& nb_mass_location = sphere_rhs->GetPosition(); real_t nb_diameter = sphere_rhs->GetDiameter(); + auto* sim = Simulation::GetActive(); + auto* param = sim->GetParam(); auto c1 = ref_mass_location; real_t r1 = 0.5 * ref_diameter; @@ -71,10 +73,10 @@ void InteractionForce::ForceBetweenSpheres(const Agent* sphere_lhs, real_t r2 = 0.5 * nb_diameter; // We take virtual bigger radii to have a distant interaction, to get a // desired density. We assume an interaction distance 3% larger then the - // agent's radius - real_t additional_radius_scaling = 1.03; - r1 *= additional_radius_scaling; - r2 *= additional_radius_scaling; + // agent's radius by default + + r1 *= param->virtual_agent_radius_scale_factor; + r2 *= param->virtual_agent_radius_scale_factor; // the 3 components of the vector c2 -> c1 real_t comp1 = c1[0] - c2[0]; real_t comp2 = c1[1] - c2[1]; diff --git a/src/core/param/param.cc b/src/core/param/param.cc index bdcecb7dc..8e48b737f 100644 --- a/src/core/param/param.cc +++ b/src/core/param/param.cc @@ -242,6 +242,8 @@ void Param::AssignFromConfig(const std::shared_ptr& config) { BDM_ASSIGN_CONFIG_VALUE(diffusion_method, "simulation.diffusion_method"); BDM_ASSIGN_CONFIG_VALUE(calculate_gradients, "simulation.calculate_gradients"); + BDM_ASSIGN_CONFIG_VALUE(virtual_agent_radius_scale_factor, + "simulation.virtual_agent_radius_scale_factor"); AssignBoundSpaceMode(config, this); AssignThreadSafetyMechanism(config, this); diff --git a/src/core/param/param.h b/src/core/param/param.h index f1d394ede..0229b7e2d 100644 --- a/src/core/param/param.h +++ b/src/core/param/param.h @@ -264,6 +264,14 @@ struct Param { /// calculate_gradients = true bool calculate_gradients = true; + /// Scale agent's virtual radius to enable distant interactions.\n + /// TOML config file: + /// Default value: `"1.03"`\n TOML + /// + /// [simulation] + /// virtual_agent_radius_scale_factor = 1.03 + real_t virtual_agent_radius_scale_factor = 1.03; + /// List of thread-safety mechanisms \n /// `kNone`: \n /// `kUserSpecified`: The user has to define all agent that must From a02895765515eddb3511de7068d9aa89907ed9fa Mon Sep 17 00:00:00 2001 From: nicogno Date: Thu, 8 Aug 2024 09:59:06 -0400 Subject: [PATCH 5/5] Fixed issue raised by SonarCloud --- src/core/interaction_force.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/interaction_force.cc b/src/core/interaction_force.cc index 27418e44d..6f03ca4d5 100644 --- a/src/core/interaction_force.cc +++ b/src/core/interaction_force.cc @@ -64,7 +64,7 @@ void InteractionForce::ForceBetweenSpheres(const Agent* sphere_lhs, real_t ref_diameter = sphere_lhs->GetDiameter(); const Real3& nb_mass_location = sphere_rhs->GetPosition(); real_t nb_diameter = sphere_rhs->GetDiameter(); - auto* sim = Simulation::GetActive(); + const auto* sim = Simulation::GetActive(); auto* param = sim->GetParam(); auto c1 = ref_mass_location;