diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a6a1db3c..38152e0e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -46,7 +46,7 @@ jobs: - name: Configure shell: bash -el {0} run: | - cmake . -Bbuild -DCMAKE_EXPORT_COMPILE_COMMANDS=ON -DCMAKE_BUILD_TYPE=Debug -DBUILD_DEMOS=ON + cmake . -Bbuild -DCMAKE_EXPORT_COMPILE_COMMANDS=ON -DCMAKE_BUILD_TYPE=Debug -DBUILD_EXAMPLES=ON -DBUILD_TESTS=ON - name: Run cppcheck on scopi shell: bash -el {0} @@ -160,6 +160,7 @@ jobs: -Bbuild \ -GNinja \ -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_EXAMPLES=ON \ -DBUILD_TESTS=ON - name: Build @@ -205,6 +206,7 @@ jobs: -Bbuild \ -GNinja \ -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_EXAMPLES=ON \ -DBUILD_TESTS=ON - name: Build diff --git a/CMakeLists.txt b/CMakeLists.txt index e8f62213..30018253 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,7 +28,6 @@ OPTION(SCOPI_USE_TBB "enable TBB" OFF) OPTION(SCOPI_USE_OPENMP "enable OpenMP" OFF) OPTION(SCOPI_USE_MOSEK "enable Mosek" OFF) OPTION(SCOPI_USE_SCS "enable SCS" OFF) -OPTION(SCOPI_USE_MKL "enable MKL" OFF) OPTION(BUILD_EXAMPLES "scopi examples" OFF) OPTION(BUILD_TESTS "scopi test suite" OFF) @@ -51,7 +50,7 @@ find_package(fmt REQUIRED) find_package(CLI11 REQUIRED) find_package(plog REQUIRED) -# find_package(nanoflann REQUIRED) +find_package(nanoflann REQUIRED) # section needed to use xtensor-blas # see https://xtensor-blas.readthedocs.io/en/latest/performance.html @@ -77,14 +76,6 @@ if(SCOPI_USE_OPENMP) find_package(OpenMP REQUIRED) endif() -if(SCOPI_USE_MKL) - if(SCOPI_USE_TBB) - set(MKL_THREAD_LAYER "TBB") - endif() - - find_package(MKL REQUIRED) -endif() - if(SCOPI_USE_MOSEK) find_package(MOSEK REQUIRED) @@ -102,16 +93,6 @@ endif() set(CMAKE_NO_SYSTEM_FROM_IMPORTED TRUE) set(SCOPI_SRC - # src/contact/contact_brute_force.cpp - # src/contact/contact_kdtree.cpp - # src/problems/DryWithoutFriction.cpp - # src/problems/DryWithFriction.cpp - # src/problems/DryWithFrictionBase.cpp - # src/problems/DryWithFrictionFixedPoint.cpp - # src/solvers/ConstraintMosek.cpp - # src/solvers/OptimUzawaBase.cpp - # src/solvers/OptimProjectedGradient.cpp - # src/solvers/projection.cpp src/vap/vap_fixed.cpp src/vap/vap_fpd.cpp src/vap/vap_projection.cpp @@ -133,7 +114,7 @@ target_link_libraries(scopi PUBLIC xtensor CLI11::CLI11 plog::plog - # nanoflann::nanoflann + nanoflann::nanoflann ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES}) # set_target_properties(scopi PROPERTIES @@ -143,10 +124,6 @@ target_link_libraries(scopi PUBLIC xtensor # VERSION ${${PROJECT_NAME}_VERSION} # SOVERSION ${VERSION_MAJOR} # OUTPUT_NAME "libscopi") -if(SCOPI_USE_MKL) - target_link_libraries(scopi PUBLIC MKL::MKL) - target_compile_definitions(scopi PUBLIC SCOPI_USE_MKL) -endif() if(SCOPI_USE_TBB) target_link_libraries(scopi PUBLIC TBB::tbb) diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index ad6119d7..c26b184c 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -1,27 +1,12 @@ set(SCOPI_EXAMPLES - # two_ellipsoids.cpp - # two_ellipsoids_asymmetrical_2d.cpp - # two_ellipsoids_asymmetrical_3d.cpp - # two_ellipsoids_spheres_2d.cpp - # two_ellipsoids_spheres_3d.cpp - # test_quentin.cpp + two_ellipsoids.cpp two_spheres.cpp - # two_spheres_periodic.cpp - # sphere_plan.cpp - # viscosity.cpp - # two_spheres_viscosity.cpp - # two_worms.cpp - # three_worms_periodic.cpp - # convergence_scheme.cpp - # box_spheres_3d.cpp - # pile_of_sand_spheres.cpp - # pile_of_sand_superellipsoids.cpp + sphere_plane.cpp + two_worms.cpp + box_spheres_3d.cpp + pile_of_sand_spheres.cpp segment.cpp - # worms_simulation.cpp ) include(generator) generate_executable(${SCOPI_EXAMPLES}) - -# add_subdirectory(critical) -# add_subdirectory(proceeding) diff --git a/demos/box_spheres_3d.cpp b/demos/box_spheres_3d.cpp index 51a12efa..886805be 100644 --- a/demos/box_spheres_3d.cpp +++ b/demos/box_spheres_3d.cpp @@ -1,29 +1,22 @@ -#include -#include +#include +#include #include -#include #include +#include #include #include -#include -#include -#include -#include -#include - int main() { -#ifdef SCOPI_USE_MKL - plog::init(plog::info, "box_spheres_3d_very_large_apgd_ar.log"); - constexpr std::size_t dim = 3; double PI = xt::numeric_constants::PI; + scopi::initialize("Spheres into a box simulation"); - std::size_t total_it = 100; + std::size_t total_it = 1000; double width_box = 10.; - std::size_t n = 50; // n^3 spheres - double g = 1.; + // std::size_t n = 50; // n^3 spheres + std::size_t n = 10; // n^3 spheres + double g = 1.; double r0 = width_box / n / 2.; double dt = 0.2 * 0.9 * r0 / (2. * g); ////// dt = 0.2*r/(sqrt(2*width_box*g)) @@ -38,31 +31,31 @@ int main() const xt::xtensor_fixed> axes_y({0., 1., 0.}); const xt::xtensor_fixed> axes_z({0., 0., 1.}); - scopi::plan p_left( + scopi::plane p_left( { {0., 0., 0.} }, {scopi::quaternion(0., axes_z)}); particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( + scopi::plane p_right( { {width_box + 2 * r0, 0., 0.} }, {scopi::quaternion(0., axes_z)}); particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( + scopi::plane p_horizontal( { {0., 0., 0.} }, {scopi::quaternion(PI / 2., axes_z)}); particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( + scopi::plane p_front( { {0., 0., 0.} }, {scopi::quaternion(PI / 2., axes_y)}); particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( + scopi::plane p_back( { {0., 0., width_box + 2 * r0} }, @@ -95,16 +88,15 @@ int main() } } - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = 99; - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; + scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles); + auto params = solver.get_params(); + // params.solver_params.output_frequency = 99; + params.optim_params.tolerance = 1e-3; + params.optim_params.alpha = rho; + params.contact_method_params.dmax = 0.9 * r0; + params.contact_method_params.kd_tree_radius = params.contact_method_params.dmax + 2. * 0.9 * r0; + + solver.run(dt, total_it); - solver.run(total_it); -#endif return 0; } diff --git a/demos/convergence_scheme.cpp b/demos/convergence_scheme.cpp deleted file mode 100644 index 4260b131..00000000 --- a/demos/convergence_scheme.cpp +++ /dev/null @@ -1,84 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - plog::init(plog::warning, "convergence_scheme.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - std::vector dt({0.1, 0.05, 0.01, 0.005, 0.001}); - std::vector total_it({100, 200, 1000, 2000, 10000}); - std::vector mu_vec({0., 0.1, 0.5, 1.}); - std::vector alpha_vec({PI / 6., PI / 4., PI / 3.}); - - for (auto mu : mu_vec) - { - for (auto alpha : alpha_vec) - { - for (std::size_t i = 0; i < dt.size(); ++i) - { - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, - dt[i]); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = mu; - solver.run(total_it[i]); - - auto pos = particles.pos(); - auto q = particles.q(); - auto tmp = scopi::analytical_solution_sphere_plan(alpha, mu, dt[i] * (total_it[i] + 1), radius, g, h); - auto pos_analytical = tmp.first; - auto q_analytical = scopi::quaternion(tmp.second); - double error_pos = xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical); - double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); - - auto v = particles.v(); - auto omega = particles.omega(); - tmp = scopi::analytical_solution_sphere_plan_velocity(alpha, mu, dt[i] * (total_it[i] + 1), radius, g, h); - auto v_analytical = tmp.first; - auto omega_analytical = tmp.second; - double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); - double error_omega = std::abs((omega(1) - omega_analytical) / omega_analytical); - - PLOG_WARNING << "mu = " << mu << " alpha = " << alpha << " dt = " << dt[i] << '\t' << error_pos << '\t' << error_q - << '\t' << error_v << '\t' << error_omega; - } - PLOG_WARNING << std::endl; - } - } - - return 0; -} diff --git a/demos/critical/2d_case.cpp b/demos/critical/2d_case.cpp index f8230229..b5ba523f 100644 --- a/demos/critical/2d_case.cpp +++ b/demos/critical/2d_case.cpp @@ -116,8 +116,8 @@ int main() // // } - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + solver.run(dt, total_it); return 0; } diff --git a/demos/critical/2d_case_no_overlap.cpp b/demos/critical/2d_case_no_overlap.cpp index bc87f804..099a6f44 100644 --- a/demos/critical/2d_case_no_overlap.cpp +++ b/demos/critical/2d_case_no_overlap.cpp @@ -72,6 +72,6 @@ int main() } } - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + solver.run(dt, total_it); } diff --git a/demos/critical/2d_case_no_velocity.cpp b/demos/critical/2d_case_no_velocity.cpp index 7d8a9ba3..c2b1e684 100644 --- a/demos/critical/2d_case_no_velocity.cpp +++ b/demos/critical/2d_case_no_velocity.cpp @@ -57,6 +57,6 @@ int main() } } - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + solver.run(dt, total_it); } diff --git a/demos/critical/2d_case_spheres.cpp b/demos/critical/2d_case_spheres.cpp index 1dd15537..9cc9ed39 100644 --- a/demos/critical/2d_case_spheres.cpp +++ b/demos/critical/2d_case_spheres.cpp @@ -59,6 +59,6 @@ int main() } } - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + solver.run(dt, total_it); } diff --git a/demos/critical/3d_case.cpp b/demos/critical/3d_case.cpp index b47dbbeb..eb76c467 100644 --- a/demos/critical/3d_case.cpp +++ b/demos/critical/3d_case.cpp @@ -95,8 +95,8 @@ int main() particles.push_back(s1, prop1); particles.push_back(s2, prop2); - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + solver.run(dt, total_it); return 0; } diff --git a/demos/critical/worms.cpp b/demos/critical/worms.cpp index 4854097d..9d10ddf6 100644 --- a/demos/critical/worms.cpp +++ b/demos/critical/worms.cpp @@ -80,8 +80,8 @@ int main() })); } - scopi::ScopiSolver> solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver> solver(particles); + solver.run(dt, total_it); return 0; } diff --git a/demos/pile_of_sand_spheres.cpp b/demos/pile_of_sand_spheres.cpp index 0f24871d..0aff0dbe 100644 --- a/demos/pile_of_sand_spheres.cpp +++ b/demos/pile_of_sand_spheres.cpp @@ -1,18 +1,14 @@ -#include -#include +#include #include -#include #include -#include +#include #include -#include -#include -#include +#include -int main() +int main(int argc, char** argv) { - plog::init(plog::info, "pile_of_sand_spheres.log"); + scopi::initialize("Pile of sand"); constexpr std::size_t dim = 3; double PI = xt::numeric_constants::PI; @@ -28,7 +24,7 @@ int main() {0., -g, 0.} }); - scopi::plan p_horizontal( + scopi::plane p_horizontal( { {0., 0., 0.} }, @@ -72,11 +68,13 @@ int main() } } - scopi::ScopiSolver, scopi::contact_brute_force, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - solver.run(total_it); + scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( + particles); + auto params = solver.get_params(); + params.default_contact_property.mu = 0.1; + + SCOPI_PARSE(argc, argv); + solver.run(dt, total_it); return 0; } diff --git a/demos/pile_of_sand_superellipsoids.cpp b/demos/pile_of_sand_superellipsoids.cpp deleted file mode 100644 index e3c06f09..00000000 --- a/demos/pile_of_sand_superellipsoids.cpp +++ /dev/null @@ -1,120 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -template -void add_obstacle(scopi::scopi_container& particles, double x, double r) -{ - // scopi::sphere s({{i*2.*r, -r}}, r); - scopi::superellipsoid s( - { - {x, -r} - }, - {r, r}, - 1.); - particles.push_back(s, scopi::property().deactivate()); -} - -int main() -{ -#ifdef SCOPI_USE_MKL - - plog::init(plog::info, "pile_of_sand_superellipsoids.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double dt = 0.01; - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 10; // n^2 ellipses - double g = 1.; - - double r = width_box / 2. / (n + 1); - double r_obs = r / 10.; - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g} - }); - - // obstacles - // scopi::plan p({{0., 0}}, PI/2.); - // particles.push_back(p, scopi::property().deactivate()); - double dist_obs = -width_box; - while (dist_obs < 2. * width_box) - { - add_obstacle(particles, dist_obs, r_obs); - dist_obs += 2 * r_obs; - } - PLOG_INFO << particles.size() << " obstacles"; - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_rx(0.5 * r, 1.5 * r); - - for (std::size_t j = 1; j < n; ++j) - { - double x = 0.; - while (x < width_box) - // for (std::size_t i = 0; i < n; ++i) - { - double m = distrib_m(generator); - double rx = distrib_rx(generator); - x += rx; - // scopi::sphere s({{x, r + j*2.*r}}, rx); - scopi::superellipsoid s( - { - {x, r + j * 2. * r} - }, - {rx, r}, - 1.); - x += rx; - particles.push_back(s, prop.mass(m).moment_inertia(m * PI / 4. * 2. * rx * r * r * r)); - } - } - - // j = 0 - double dec_x = 0.5 * r; - double x = 0.; - // for (std::size_t i = 0; i < n; ++i) - while (x < width_box) - { - double m = distrib_m(generator); - double rx = distrib_rx(generator); - x += rx; - // scopi::sphere s({{x + dec_x, r}}, rx); - scopi::superellipsoid s( - { - {x + dec_x, r} - }, - {rx, r}, - 1.); - x += rx; - particles.push_back(s, prop.mass(m).moment_inertia(m * PI / 4. * 2. * rx * r * r * r)); - } - - scopi::ScopiSolver, scopi::contact_brute_force, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-6; - params.contact_params.dmax = 2. * 1.5 * r; - // params.optim_params.rho = 2.; - // params.solver_params.output_frequency = 2; - // params.optim_params.change_default_tol_mosek = false; - // params.problem_params.mu = 0.1; - solver.run(total_it); -#endif - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/CMakeLists.txt b/demos/proceeding/disk_on_inclined_plane/CMakeLists.txt deleted file mode 100644 index 842b74c3..00000000 --- a/demos/proceeding/disk_on_inclined_plane/CMakeLists.txt +++ /dev/null @@ -1,29 +0,0 @@ -list(APPEND SCOPI_PROCEEDING_DISK_ON_INCLINED_PLANE - distance_and_velocity_with_friction_mu01_convex.cpp - distance_and_velocity_with_friction_mu01_fixedPoint.cpp - distance_and_velocity_with_friction_mu1_convex.cpp - distance_and_velocity_with_friction_mu1_fixedPoint.cpp - error_with_friction_mu01_convex.cpp - error_with_friction_mu01_fixedPoint.cpp - error_with_friction_mu1_convex.cpp - error_with_friction_mu1_fixedPoint.cpp -) - -if(SCOPI_USE_MKL) - list(APPEND SCOPI_PROCEEDING_DISK_ON_INCLINED_PLANE - distance_and_velocity_without_friction.cpp - oneStep_pgd_smallTol.cpp - oneStep_apgd_smallTol.cpp - oneStep_apgdas_smallTol.cpp - oneStep_apgdar_smallTol.cpp - oneStep_apgdasr_smallTol.cpp - oneStep_pgd_largeTol.cpp - oneStep_apgd_largeTol.cpp - oneStep_apgdas_largeTol.cpp - oneStep_apgdar_largeTol.cpp - oneStep_apgdasr_largeTol.cpp - error_without_friction.cpp - ) -endif() - -generate_executable(${SCOPI_PROCEEDING_DISK_ON_INCLINED_PLANE}) diff --git a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu01_convex.cpp b/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu01_convex.cpp deleted file mode 100644 index 622ee9dc..00000000 --- a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu01_convex.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 5: Disk falling on a plane with friction. Distance, tangential velocity, normal velocity, and rotation. - // mu = 0.1, convex scheme - plog::init(plog::error, "distance_and_velocity_with_friction_mu01_convex.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 200; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.problem_params.mu = 0.1; - params.optim_params.change_default_tol_mosek = false; - params.solver_params.write_velocity = true; - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu01_fixedPoint.cpp b/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu01_fixedPoint.cpp deleted file mode 100644 index ed9c9e55..00000000 --- a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu01_fixedPoint.cpp +++ /dev/null @@ -1,58 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 5: Disk falling on a plane with friction. Distance, tangential velocity, normal velocity, and rotation. - // mu = 0.1, fixed point algorithm - plog::init(plog::error, "distance_and_velocity_with_friction_mu01_convex.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 200; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.problem_params.mu = 0.1; - params.problem_params.tol_fixed_point = 1e-6; - params.optim_params.change_default_tol_mosek = false; - params.solver_params.write_velocity = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu1_convex.cpp b/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu1_convex.cpp deleted file mode 100644 index ed263683..00000000 --- a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu1_convex.cpp +++ /dev/null @@ -1,57 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 5: Disk falling on a plane with friction. Distance, tangential velocity, normal velocity, and rotation. - // mu = 1, convex scheme - plog::init(plog::error, "distance_and_velocity_with_friction_mu01_convex.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 200; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.problem_params.mu = 1.; - params.optim_params.change_default_tol_mosek = false; - params.solver_params.write_velocity = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu1_fixedPoint.cpp b/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu1_fixedPoint.cpp deleted file mode 100644 index b269d506..00000000 --- a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_with_friction_mu1_fixedPoint.cpp +++ /dev/null @@ -1,58 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 5: Disk falling on a plane with friction. Distance, tangential velocity, normal velocity, and rotation. - // mu = 1, fixed point algorithm - plog::init(plog::error, "distance_and_velocity_with_friction_mu01_convex.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 200; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.problem_params.mu = 1.; - params.problem_params.tol_fixed_point = 1e-6; - params.optim_params.change_default_tol_mosek = false; - params.solver_params.write_velocity = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_without_friction.cpp b/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_without_friction.cpp deleted file mode 100644 index ac91f933..00000000 --- a/demos/proceeding/disk_on_inclined_plane/distance_and_velocity_without_friction.cpp +++ /dev/null @@ -1,57 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 3: Disk falling on a plane without friction. Distance, tangential velocity, and normal velocity. - plog::init(plog::error, "distance_and_velocity_without_friction.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 200; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.solver_params.write_velocity = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu01_convex.cpp b/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu01_convex.cpp deleted file mode 100644 index 4e03302a..00000000 --- a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu01_convex.cpp +++ /dev/null @@ -1,80 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 6: Disk falling on a plane with friction. L^2 error on the position of the center and the rotation of the disk as a function - // of dt. mu = 0.1, convex scheme - plog::init(plog::warning, "error_with_friction_mu01_convex.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - double alpha = PI / 6.; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - std::vector dt({0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05}); - std::vector total_it({100000, 50000, 20000, 10000, 5000, 2000, 1000, 500, 200}); - - for (std::size_t i = 0; i < dt.size(); ++i) - { - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - double error_pos = 0.; - double error_rot = 0.; - for (std::size_t n = 1; n < total_it[i]; ++n) - { - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, - dt[i]); - auto params = solver.get_params(); - params.problem_params.mu = 0.1; - params.solver_params.output_frequency = std::size_t(-1); - solver.run(n, n - 1); - - auto tmp = scopi::analytical_solution_sphere_plan(alpha, params.problem_params.mu, dt[i] * n, radius, g, h); - - auto pos = particles.pos(); - auto pos_analytical = tmp.first; - error_pos += (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)) - * (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)); - - auto q = particles.q(); - auto q_analytical = scopi::quaternion(tmp.second); - error_rot += (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)) - * (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)); - } - PLOG_WARNING << "dt = " << dt[i] << " err pos = " << std::sqrt(dt[i] * error_pos) << " err rot = " << std::sqrt(dt[i] * error_rot); - } - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu01_fixedPoint.cpp b/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu01_fixedPoint.cpp deleted file mode 100644 index 948e993a..00000000 --- a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu01_fixedPoint.cpp +++ /dev/null @@ -1,83 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 6: Disk falling on a plane with friction. L^2 error on the position of the center and the rotation of the disk as a function - // of dt. mu = 0.1, fixed point algorithm - plog::init(plog::warning, "error_with_friction_mu01_fixedPoint.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - double alpha = PI / 6.; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - std::vector dt({0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05}); - std::vector total_it({100000, 50000, 20000, 10000, 5000, 2000, 1000, 500, 200}); - - for (std::size_t i = 0; i < dt.size(); ++i) - { - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - double error_pos = 0.; - double error_rot = 0.; - for (std::size_t n = 1; n < total_it[i]; ++n) - { - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt[i]); - auto params = solver.get_params(); - params.problem_params.mu = 0.1; - params.problem_params.tol_fixed_point = 1e-2; - params.solver_params.output_frequency = -1; - - solver.run(n, n - 1); - - auto tmp = scopi::analytical_solution_sphere_plan(alpha, params.problem_params.mu, dt[i] * n, radius, g, h); - - auto pos = particles.pos(); - auto pos_analytical = tmp.first; - error_pos += (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)) - * (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)); - - auto q = particles.q(); - auto q_analytical = scopi::quaternion(tmp.second); - error_rot += (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)) - * (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)); - } - PLOG_WARNING << "dt = " << dt[i] << " err pos = " << std::sqrt(dt[i] * error_pos) << " err rot = " << std::sqrt(dt[i] * error_rot); - } - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu1_convex.cpp b/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu1_convex.cpp deleted file mode 100644 index 3d837755..00000000 --- a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu1_convex.cpp +++ /dev/null @@ -1,81 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 6: Disk falling on a plane with friction. L^2 error on the position of the center and the rotation of the disk as a function - // of dt. mu = 1, convex scheme - plog::init(plog::warning, "error_with_friction_mu1_convex.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - double alpha = PI / 6.; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - std::vector dt({0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05}); - std::vector total_it({100000, 50000, 20000, 10000, 5000, 2000, 1000, 500, 200}); - - for (std::size_t i = 0; i < dt.size(); ++i) - { - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - double error_pos = 0.; - double error_rot = 0.; - for (std::size_t n = 1; n < total_it[i]; ++n) - { - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, - dt[i]); - auto params = solver.get_params(); - params.problem_params.mu = 1.; - params.solver_params.output_frequency = -1; - - solver.run(n, n - 1); - - auto tmp = scopi::analytical_solution_sphere_plan(alpha, params.problem_params.mu, dt[i] * n, radius, g, h); - - auto pos = particles.pos(); - auto pos_analytical = tmp.first; - error_pos += (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)) - * (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)); - - auto q = particles.q(); - auto q_analytical = scopi::quaternion(tmp.second); - error_rot += (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)) - * (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)); - } - PLOG_WARNING << "dt = " << dt[i] << " err pos = " << std::sqrt(dt[i] * error_pos) << " err rot = " << std::sqrt(dt[i] * error_rot); - } - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu1_fixedPoint.cpp b/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu1_fixedPoint.cpp deleted file mode 100644 index d44260ef..00000000 --- a/demos/proceeding/disk_on_inclined_plane/error_with_friction_mu1_fixedPoint.cpp +++ /dev/null @@ -1,83 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 6: Disk falling on a plane with friction. L^2 error on the position of the center and the rotation of the disk as a function - // of dt. mu = 1, fixed point algorithm - plog::init(plog::warning, "error_with_friction_mu1_fixedPoint.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - double alpha = PI / 6.; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - std::vector dt({0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05}); - std::vector total_it({100000, 50000, 20000, 10000, 5000, 2000, 1000, 500, 200}); - - for (std::size_t i = 0; i < dt.size(); ++i) - { - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - double error_pos = 0.; - double error_rot = 0.; - for (std::size_t n = 1; n < total_it[i]; ++n) - { - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt[i]); - auto params = solver.get_params(); - params.problem_params.mu = 1.; - params.problem_params.tol_fixed_point = 1e-2; - params.solver_params.output_frequency = -1; - - solver.run(n, n - 1); - - auto tmp = scopi::analytical_solution_sphere_plan(alpha, params.problem_params.mu, dt[i] * n, radius, g, h); - - auto pos = particles.pos(); - auto pos_analytical = tmp.first; - error_pos += (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)) - * (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)); - - auto q = particles.q(); - auto q_analytical = scopi::quaternion(tmp.second); - error_rot += (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)) - * (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)); - } - PLOG_WARNING << "dt = " << dt[i] << " err pos = " << std::sqrt(dt[i] * error_pos) << " err rot = " << std::sqrt(dt[i] * error_rot); - } - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/error_without_friction.cpp b/demos/proceeding/disk_on_inclined_plane/error_without_friction.cpp deleted file mode 100644 index d08e2931..00000000 --- a/demos/proceeding/disk_on_inclined_plane/error_without_friction.cpp +++ /dev/null @@ -1,82 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -int main() -{ - // Figure 4: Disk falling on a plane without friction. L^2 error on the position of the center of the disk as a function of dt. - plog::init(plog::warning, "error_without_friction.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = 2. * radius; - double alpha = PI / 6.; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - std::vector dt({0.0001, 0.0002, 0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05}); - std::vector total_it({100000, 50000, 20000, 10000, 5000, 2000, 1000, 500, 200}); - - for (std::size_t i = 0; i < dt.size(); ++i) - { - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - double error_pos = 0.; - double error_rot = 0.; - for (std::size_t n = 1; n < total_it[i]; ++n) - { - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt[i]); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.solver_params.output_frequency = 100000; - - solver.run(n, n - 1); - - auto tmp = scopi::analytical_solution_sphere_plan(alpha, 0., dt[i] * n, radius, g, h); - - auto pos = particles.pos(); - auto pos_analytical = tmp.first; - error_pos += (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)) - * (xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical)); - - auto q = particles.q(); - auto q_analytical = scopi::quaternion(tmp.second); - error_rot += (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)) - * (xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical)); - } - PLOG_WARNING << "dt = " << dt[i] << " err pos = " << std::sqrt(dt[i] * error_pos) << " err rot = " << std::sqrt(dt[i] * error_rot); - } - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgd_largeTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgd_largeTol.cpp deleted file mode 100644 index e2d7f1d8..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgd_largeTol.cpp +++ /dev/null @@ -1,57 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 1: Disk placed on a plane without friction. Number of iterations for the APGD algorithm and tolerance 10^{-3}. - plog::init(plog::info, "disk_on_inclined_plane_oneStep_apdg_largeTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = 2.; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgd_smallTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgd_smallTol.cpp deleted file mode 100644 index f4f9d148..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgd_smallTol.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 2: Disk placed on a plane without friction. Evolution of the constraint and cost as a function of the iterations for the APGD - // algorithm. Table 1: Disk placed on a plane without friction. Number of iterations for the APGD algorithm and tolerance 10^{-9}. The - // constraint is the first column in the log, the cost is the second column. - plog::init(plog::verbose, "disk_on_inclined_plane_oneStep_apdg_smallTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.optim_params.verbose = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdar_largeTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgdar_largeTol.cpp deleted file mode 100644 index 24148dd1..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdar_largeTol.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 1: Disk placed on a plane without friction. Number of iterations for the APGD-AR algorithm and tolerance 10^{-3}. - plog::init(plog::info, "disk_on_inclined_plane_oneStep_apdgar_largeTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = 2.; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdar_smallTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgdar_smallTol.cpp deleted file mode 100644 index ee9225bf..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdar_smallTol.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 2: Disk placed on a plane without friction. Evolution of the constraint and cost as a function of the iterations for the - // APGD-AR algorithm. Table 1: Disk placed on a plane without friction. Number of iterations for the APGD-AR algorithm and tolerance - // 10^{-9}. The constraint is the first column in the log, the cost is the second column. - plog::init(plog::verbose, "disk_on_inclined_plane_oneStep_apdgar_smallTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.optim_params.verbose = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdas_largeTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgdas_largeTol.cpp deleted file mode 100644 index 7625f34b..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdas_largeTol.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 1: Disk placed on a plane without friction. Number of iterations for the APGD-AS algorithm and tolerance 10^{-3}. - plog::init(plog::info, "disk_on_inclined_plane_oneStep_apdgas_largeTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = 2.; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdas_smallTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgdas_smallTol.cpp deleted file mode 100644 index 60c330e8..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdas_smallTol.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 2: Disk placed on a plane without friction. Evolution of the constraint and cost as a function of the iterations for the - // APGD-AS algorithm. Table 1: Disk placed on a plane without friction. Number of iterations for the APGD-AS algorithm and tolerance - // 10^{-9}. The constraint is the first column in the log, the cost is the second column. - plog::init(plog::verbose, "disk_on_inclined_plane_oneStep_apdgas_smallTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.optim_params.verbose = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdasr_largeTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgdasr_largeTol.cpp deleted file mode 100644 index 2a3c886a..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdasr_largeTol.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 1: Disk placed on a plane without friction. Number of iterations for the APGD-ASR algorithm and tolerance 10^{-3}. - plog::init(plog::info, "disk_on_inclined_plane_oneStep_apdgasr_largeTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = 2.; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdasr_smallTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_apgdasr_smallTol.cpp deleted file mode 100644 index 7de74a8e..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_apgdasr_smallTol.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 2: Disk placed on a plane without friction. Evolution of the constraint and cost as a function of the iterations for the - // APGD-ASR algorithm. Table 1: Disk placed on a plane without friction. Number of iterations for the APGD-ASR algorithm and tolerance - // 10^{-9}. The constraint is the first column in the log, the cost is the second column. - plog::init(plog::verbose, "disk_on_inclined_plane_oneStep_apdgasr_smallTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.optim_params.verbose = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_pgd_largeTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_pgd_largeTol.cpp deleted file mode 100644 index 5953813d..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_pgd_largeTol.cpp +++ /dev/null @@ -1,57 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 1: Disk placed on a plane without friction. Number of iterations for the PGD algorithm and tolerance 10^{-3}. - plog::init(plog::info, "disk_on_inclined_plane_oneStep_pdg_largeTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = 2.; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/disk_on_inclined_plane/oneStep_pgd_smallTol.cpp b/demos/proceeding/disk_on_inclined_plane/oneStep_pgd_smallTol.cpp deleted file mode 100644 index 9b6e9038..00000000 --- a/demos/proceeding/disk_on_inclined_plane/oneStep_pgd_smallTol.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 2: Disk placed on a plane without friction. Evolution of the constraint and cost as a function of the iterations for the PGD - // algorithm. Table 1: Disk placed on a plane without friction. Number of iterations for the PGD algorithm and tolerance 10^{-9}. The - // constraint is the first column in the log, the cost is the second column. - plog::init(plog::verbose, "disk_on_inclined_plane_oneStep_pdg_smallTol.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double mass = 1.; - double h = radius; - auto prop = scopi::property().mass(mass).moment_inertia(mass * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 1; - double alpha = PI / 6.; - - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2. - alpha); - scopi::sphere s( - { - {h * std::sin(alpha), h * std::cos(alpha)} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {0., -g} - })); - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-9; - params.optim_params.rho = 2.; - params.optim_params.verbose = true; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/CMakeLists.txt b/demos/proceeding/multiparticles_tests/CMakeLists.txt deleted file mode 100644 index 100ae9b6..00000000 --- a/demos/proceeding/multiparticles_tests/CMakeLists.txt +++ /dev/null @@ -1,32 +0,0 @@ -list(APPEND SCOPI_PROCEEDING_MULTIPARTICLES - pile_of_sand_friction_small_config_mu01_convex.cpp - pile_of_sand_friction_small_config_mu01_fixed_point.cpp - pile_of_sand_friction_small_config_mu1_convex.cpp - pile_of_sand_friction_small_config_mu1_fixed_point.cpp - pile_of_sand_friction_large_config_mu01_convex.cpp - pile_of_sand_friction_large_config_mu01_fixed_point.cpp - pile_of_sand_friction_large_config_mu1_convex.cpp - pile_of_sand_friction_large_config_mu1_fixed_point.cpp - pile_of_sand_friction_figure_mu01_convex.cpp - pile_of_sand_friction_figure_mu01_fixed_point.cpp - pile_of_sand_friction_figure_mu1_convex.cpp - pile_of_sand_friction_figure_mu1_fixed_point.cpp -) - -if(SCOPI_USE_MKL) - list(APPEND SCOPI_PROCEEDING_MULTIPARTICLES - spheres_in_box_small_config_pgd.cpp - spheres_in_box_small_config_apgd.cpp - spheres_in_box_small_config_apgdas.cpp - spheres_in_box_small_config_apgdar.cpp - spheres_in_box_small_config_apgdasr.cpp - spheres_in_box_large_config_pgd.cpp - spheres_in_box_large_config_apgd.cpp - spheres_in_box_large_config_apgdas.cpp - spheres_in_box_large_config_apgdar.cpp - spheres_in_box_large_config_apgdasr.cpp - pile_of_sand_ellipses.cpp - ) -endif() - -generate_executable(${SCOPI_PROCEEDING_MULTIPARTICLES}) diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_ellipses.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_ellipses.cpp deleted file mode 100644 index be82149b..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_ellipses.cpp +++ /dev/null @@ -1,108 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -template -void add_obstacle(scopi::scopi_container& particles, double x, double r) -{ - scopi::superellipsoid s( - { - {x, -r} - }, - {r, r}, - 1.); - particles.push_back(s, scopi::property().deactivate()); -} - -int main() -{ - // Figure 10: ellipses falling on a plane. - plog::init(plog::info, "pile_of_sand_ellipses.log"); - - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 10; // n^2 ellipses - double g = 5.; - - double r = width_box / 2. / (n + 1); - double r_obs = r / 10.; - double dt = 0.2 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g} - }); - - // obstacles - double dist_obs = -width_box; - while (dist_obs < 2. * width_box) - { - add_obstacle(particles, dist_obs, r_obs); - dist_obs += 2 * r_obs; - } - PLOG_INFO << particles.size() << " obstacles"; - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_rx(0.5 * r, 1.5 * r); - - for (std::size_t j = 1; j < n; ++j) - { - double x = 0.; - while (x < width_box) - { - double m = distrib_m(generator); - double rx = distrib_rx(generator); - x += rx; - scopi::superellipsoid s( - { - {x, r + j * 2. * r} - }, - {rx, r}, - 1.); - x += rx; - particles.push_back(s, prop.mass(m).moment_inertia(m * PI / 4. * 2. * rx * r * r * r)); - } - } - - // j = 0 - double dec_x = 0.5 * r; - double x = 0.; - while (x < width_box) - { - double m = distrib_m(generator); - double rx = distrib_rx(generator); - x += rx; - scopi::superellipsoid s( - { - {x + dec_x, r} - }, - {rx, r}, - 1.); - x += rx; - particles.push_back(s, prop.mass(m).moment_inertia(m * PI / 4. * 2. * rx * r * r * r)); - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = 0.2 / dt / dt; - params.solver_params.path = "pile_sand_friction"; - params.solver_params.filename = "ellipses/scopi_objects_"; - params.contact_params.dmax = 2. * r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 3. * r; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu01_convex.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu01_convex.cpp deleted file mode 100644 index d8d77355..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu01_convex.cpp +++ /dev/null @@ -1,90 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 9: spheres falling on a plane with friction. - // mu = 0.1, convex scheme. - plog::init(plog::info, "pile_of_sand_spheres_figure_mu01_convex.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 10; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - // params.problem_params.mu = 0.1; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "nofriction_convex_"; - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu01_fixed_point.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu01_fixed_point.cpp deleted file mode 100644 index 2c147497..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu01_fixed_point.cpp +++ /dev/null @@ -1,92 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 9: spheres falling on a plane with friction. - // mu = 0.1, fixed point algorithm. - plog::init(plog::info, "pile_of_sand_spheres_figure_mu01_fixed_point.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 10; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - params.problem_params.tol_fixed_point = 1e-2; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "mu01_fixed_point_"; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu1_convex.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu1_convex.cpp deleted file mode 100644 index 05448d04..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu1_convex.cpp +++ /dev/null @@ -1,90 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 9: spheres falling on a plane with friction. - // mu = 1, convex scheme. - plog::init(plog::info, "pile_of_sand_spheres_figure_mu1_convex.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 10; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 1.; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_friction"; - params.solver_params.filename = "mu1_convex_"; - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu1_fixed_point.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu1_fixed_point.cpp deleted file mode 100644 index bf6f280c..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_figure_mu1_fixed_point.cpp +++ /dev/null @@ -1,92 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 9: spheres falling on a plane with friction. - // mu = 1, fixed point algorithm. - plog::init(plog::info, "pile_of_sand_spheres_figure_mu1_fixed_point.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 10; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 1.; - params.problem_params.tol_fixed_point = 1e-2; - params.problem_params.max_iter_fixed_point = 100; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "mu1_fixed_point_"; - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu01_convex.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu01_convex.cpp deleted file mode 100644 index 2ffb4984..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu01_convex.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 12^3 spheres falling on a plane with friction. - // mu = 0.1, convex scheme. - plog::init(plog::info, "pile_of_sand_spheres_large_config_mu01_convex.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 12; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "mu01_convex_"; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu01_fixed_point.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu01_fixed_point.cpp deleted file mode 100644 index 583fb858..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu01_fixed_point.cpp +++ /dev/null @@ -1,92 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 12^3 spheres falling on a plane with friction. - // mu = 0.1, fixed point algorithm. - plog::init(plog::info, "pile_of_sand_spheres_large_config_mu01_fixed_point.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 12; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - params.problem_params.tol_fixed_point = 1e-2; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "mu01_fixed_point_"; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu1_convex.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu1_convex.cpp deleted file mode 100644 index 26e221c4..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu1_convex.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 12^3 spheres falling on a plane with friction. - // mu = 1, convex scheme. - plog::init(plog::info, "pile_of_sand_spheres_large_config_mu1_convex.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 12; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 1.; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "mu1_convex_"; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu1_fixed_point.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu1_fixed_point.cpp deleted file mode 100644 index 39d1d7e1..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_large_config_mu1_fixed_point.cpp +++ /dev/null @@ -1,93 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 12^3 spheres falling on a plane with friction. - // mu = 1, fixed point algorithm. - plog::init(plog::info, "pile_of_sand_spheres_large_config_mu1_fixed_point.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 12; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 1.; - params.problem_params.tol_fixed_point = 1e-2; - params.problem_params.max_iter_fixed_point = 100; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = 20; - params.solver_params.path = "pile_sand_no_friction"; - params.solver_params.filename = "mu1_fixed_point_"; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu01_convex.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu01_convex.cpp deleted file mode 100644 index 42d967ce..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu01_convex.cpp +++ /dev/null @@ -1,89 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 8^3 spheres falling on a plane with friction. - // mu = 0.1, convex scheme. - plog::init(plog::info, "pile_of_sand_spheres_small_config_mu01_convex.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = std::size_t(-1); - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu01_fixed_point.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu01_fixed_point.cpp deleted file mode 100644 index 73b5f964..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu01_fixed_point.cpp +++ /dev/null @@ -1,90 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 8^3 spheres falling on a plane with friction. - // mu = 0.1, fixed point algorithm. - plog::init(plog::info, "pile_of_sand_spheres_small_config_mu01_fixed_point.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - params.problem_params.tol_fixed_point = 1e-2; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = std::size_t(-1); - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu1_convex.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu1_convex.cpp deleted file mode 100644 index 3508fe93..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu1_convex.cpp +++ /dev/null @@ -1,89 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 8^3 spheres falling on a plane with friction. - // mu = 1, convex scheme. - plog::init(plog::info, "pile_of_sand_spheres_small_config_mu1_convex.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 1.; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = std::size_t(-1); - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu1_fixed_point.cpp b/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu1_fixed_point.cpp deleted file mode 100644 index 92484f81..00000000 --- a/demos/proceeding/multiparticles_tests/pile_of_sand_friction_small_config_mu1_fixed_point.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 3: 8^3 spheres falling on a plane with friction. - // mu = 1, fixed point algorithm. - plog::init(plog::info, "pile_of_sand_spheres_small_config_mu1_fixed_point.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - std::size_t total_it = 1000; - double g = 1.; - - double r = width_box / 2. / (n + 1); - double dt = 0.1 * r / (std::sqrt(2. * width_box * g)); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - PI / 2.); - particles.push_back(p_horizontal, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_m(1., 2.); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 1; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r, r + j * 2. * r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - // j = 0 - double dec_x = 0.5 * r; - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t k = 0; k < n; ++k) - { - double m = distrib_m(generator); - scopi::sphere s( - { - {i * 2. * r + dec_x, r, k * 2. * r} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 1.; - params.problem_params.tol_fixed_point = 1e-2; - params.problem_params.max_iter_fixed_point = 100; - params.contact_params.dmax = r; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * r; - params.solver_params.output_frequency = std::size_t(-1); - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgd.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgd.cpp deleted file mode 100644 index 4ad1e1fc..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgd.cpp +++ /dev/null @@ -1,110 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 8: number of iterations vs number of active contatcs for spheres falling in a box. - // Table 2: 50^3 spheres falling in a box. - // APGD algorithm. - plog::init(plog::info, "spheres_in_box_large_config_apgd.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 50; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdar.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdar.cpp deleted file mode 100644 index 320ed613..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdar.cpp +++ /dev/null @@ -1,113 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 7: spheres falling in a box. - // Figure 8: number of iterations vs number of active contatcs for spheres falling in a box. - // Table 2: 50^3 spheres falling in a box. - // APGD-AR algorithm. - plog::init(plog::info, "spheres_in_box_large_config_apgdar.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - // std::size_t total_it = 1000; // For the number of contacts, like the other solvers - std::size_t total_it = 2000; // For output of stationary state - double width_box = 10.; - std::size_t n = 50; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = 100; - params.solver_params.path = "proceeding"; - params.solver_params.filename = "spheres_in_box"; - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdas.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdas.cpp deleted file mode 100644 index d0d3cda6..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdas.cpp +++ /dev/null @@ -1,109 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 8: number of iterations vs number of active contatcs for spheres falling in a box. - // Table 2: 50^3 spheres falling in a box. - // APGD-AS algorithm. - plog::init(plog::info, "spheres_in_box_large_config_apgdas.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 50; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdasr.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdasr.cpp deleted file mode 100644 index d92711c6..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_apgdasr.cpp +++ /dev/null @@ -1,109 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 8: number of iterations vs number of active contatcs for spheres falling in a box. - // Table 2: 50^3 spheres falling in a box. - // APGD-ASR algorithm. - plog::init(plog::info, "spheres_in_box_large_config_apgdasr.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 50; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = 2000; - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_pgd.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_pgd.cpp deleted file mode 100644 index 274c23e6..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_large_config_pgd.cpp +++ /dev/null @@ -1,110 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Figure 8: number of iterations vs number of active contatcs for spheres falling in a box. - // Table 2: 50^3 spheres falling in a box. - // PGD algorithm. - plog::init(plog::info, "spheres_in_box_large_config_pgd.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 50; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgd.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgd.cpp deleted file mode 100644 index 67733257..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgd.cpp +++ /dev/null @@ -1,109 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 2: 8^3 spheres falling in a box. - // APGD algorithm. - plog::init(plog::info, "spheres_in_box_small_config_apgd.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdar.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdar.cpp deleted file mode 100644 index 21d1d7f0..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdar.cpp +++ /dev/null @@ -1,108 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 2: 8^3 spheres falling in a box. - // APGD-AR algorithm. - plog::init(plog::info, "spheres_in_box_small_config_apgdar.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdas.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdas.cpp deleted file mode 100644 index 908a3667..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdas.cpp +++ /dev/null @@ -1,108 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 2: 8^3 spheres falling in a box. - // APGD-AS algorithm. - plog::init(plog::info, "spheres_in_box_small_config_apgdas.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdasr.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdasr.cpp deleted file mode 100644 index ad0627bb..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_apgdasr.cpp +++ /dev/null @@ -1,110 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 2: 8^3 spheres falling in a box. - // APGD-ASR algorithm. - plog::init(plog::info, "spheres_in_box_small_config_apgdasr.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> - solver(particles, dt); - - scopi::Params, scopi::contact_kdtree, scopi::vap_fpd> params; - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_pgd.cpp b/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_pgd.cpp deleted file mode 100644 index f3e632dc..00000000 --- a/demos/proceeding/multiparticles_tests/spheres_in_box_small_config_pgd.cpp +++ /dev/null @@ -1,109 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -int main() -{ - // Table 2: 8^3 spheres falling in a box. - // PGD algorithm. - plog::init(plog::info, "spheres_in_box_small_config_pgd.log"); - - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - - std::size_t total_it = 1000; - double width_box = 10.; - std::size_t n = 8; // n^3 spheres - double g = 1.; - - double r0 = width_box / n / 2.; - double dt = 0.2 * 0.9 * r0 / (std::sqrt(2. * width_box * g)); - double rho = 0.2 / (dt * dt); - - scopi::scopi_container particles; - auto prop = scopi::property().force({ - {0., -g, 0.} - }); - - const xt::xtensor_fixed> axes_y({0., 1., 0.}); - const xt::xtensor_fixed> axes_z({0., 0., 1.}); - - scopi::plan p_left( - { - {0., 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_left, scopi::property().deactivate()); - scopi::plan p_right( - { - {width_box + 2 * r0, 0., 0.} - }, - {scopi::quaternion(0., axes_z)}); - particles.push_back(p_right, scopi::property().deactivate()); - scopi::plan p_horizontal( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_z)}); - particles.push_back(p_horizontal, scopi::property().deactivate()); - scopi::plan p_front( - { - {0., 0., 0.} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_front, scopi::property().deactivate()); - scopi::plan p_back( - { - {0., 0., width_box + 2 * r0} - }, - {scopi::quaternion(PI / 2., axes_y)}); - particles.push_back(p_back, scopi::property().deactivate()); - - std::default_random_engine generator; - std::uniform_real_distribution distrib_r(0.8 * r0, 0.9 * r0); - std::uniform_real_distribution distrib_m(1., 2.); - std::uniform_real_distribution distrib_mouve_center(-0.05 * r0, 0.05 * r0); - - for (std::size_t i = 0; i < n; ++i) - { - for (std::size_t j = 0; j < n; ++j) - { - for (std::size_t k = 0; k < n; ++k) - { - double r = distrib_r(generator); - double m = distrib_m(generator); - double dx = distrib_mouve_center(generator); - double dy = distrib_mouve_center(generator); - double dz = distrib_mouve_center(generator); - scopi::sphere s( - { - {r0 + 0.1 + i * 2. * r0 + dx, 1.1 * r0 + 2.1 * r0 * j + dy, r0 + 0.1 + k * 2. * r0 + dz} - }, - r); - particles.push_back(s, prop.mass(m).moment_inertia({m * r * r / 2., m * r * r / 2., m * r * r / 2.})); - } - } - } - - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( - particles, - dt); - auto params = solver.get_params(); - params.solver_params.output_frequency = std::size_t(-1); - params.optim_params.tol_l = 1e-3; - params.optim_params.rho = rho; - params.contact_params.dmax = 0.9 * r0; - params.contact_params.kd_tree_radius = params.contact_params.dmax + 2. * 0.9 * r0; - - solver.run(total_it); - - return 0; -} diff --git a/demos/segment.cpp b/demos/segment.cpp index 769cc3eb..d42f91a1 100644 --- a/demos/segment.cpp +++ b/demos/segment.cpp @@ -1,41 +1,29 @@ #include -#include - -#include -#include #include #include -#include #include -#include - -#include -#include -#include -#include int main(int argc, char** argv) { constexpr std::size_t dim = 2; double dt = .001; const double max_radius = 0.06; - std::size_t total_it = 100; - std::size_t n_parts = 20; - - static plog::ColorConsoleAppender consoleAppender; - plog::init(plog::info, &consoleAppender); + std::size_t total_it = 1000; + std::size_t n_parts = 50; - CLI::App app("two segments"); + scopi::initialize("spheres passing between two segments"); + auto& app = scopi::get_app(); app.add_option("--nparts", n_parts, "Number of particles")->capture_default_str(); app.add_option("--nite", total_it, "Number of iterations")->capture_default_str(); app.add_option("--dt", dt, "Time step")->capture_default_str(); - CLI11_PARSE(app, argc, argv); + + scopi::scopi_container particles; + scopi::ScopiSolver solver(particles); + SCOPI_PARSE(argc, argv); scopi::segment seg1(scopi::type::position_t{0., 1.}, scopi::type::position_t{0.4, 1.}); scopi::segment seg2(scopi::type::position_t{0.6, 1.}, scopi::type::position_t{1., 1.}); - - scopi::scopi_container particles; particles.push_back(seg1, scopi::property().deactivate()); particles.push_back(seg2, scopi::property().deactivate()); @@ -64,22 +52,11 @@ int main(int argc, char** argv) prop); } - using problem_t = scopi::NoFriction; - using optim_solver = scopi::OptimGradient; - using contact_t = scopi::contact_kdtree; - using vap_t = scopi::vap_fixed; - using solver_type = scopi::ScopiSolver; - solver_type solver(particles, dt); - - auto params = solver.get_params(); - // params.optim_params.alpha = 0.2 / (dt * dt); - params.optim_params.dynamic_descent = true; - params.contact_params.dmax = 4 * max_radius; - - // solver.init_options(app); - // CLI11_PARSE(app, argc, argv); + auto params = solver.get_params(); + params.contact_method_params.dmax = 2 * dt; + params.contact_method_params.kd_tree_radius = 4 * max_radius; - solver.run(total_it); + solver.run(dt, total_it); return 0; } diff --git a/demos/sphere_plan.cpp b/demos/sphere_plane.cpp similarity index 64% rename from demos/sphere_plan.cpp rename to demos/sphere_plane.cpp index d7d57234..beb532e8 100644 --- a/demos/sphere_plan.cpp +++ b/demos/sphere_plane.cpp @@ -1,17 +1,13 @@ -#include +#include #include -#include +#include #include -#include -#include - -#include -#include #include +#include -int main() +int main(int argc, char** argv) { - plog::init(plog::info, "sphere_plan.log"); + scopi::initialize("sphere plane simulation"); constexpr std::size_t dim = 2; double PI = xt::numeric_constants::PI; @@ -27,7 +23,7 @@ int main() double alpha = PI / 6.; scopi::scopi_container particles; - scopi::plan p( + scopi::plane p( { {0., 0.} }, @@ -43,10 +39,13 @@ int main() {0., -g} })); - scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.problem_params.mu = 0.1; - solver.run(total_it); + scopi::ScopiSolver, scopi::contact_kdtree, scopi::vap_fpd> solver( + particles); + auto params = solver.get_params(); + params.default_contact_property.mu = 0.1; + + SCOPI_PARSE(argc, argv); + solver.run(dt, total_it); return 0; } diff --git a/demos/test_quentin.cpp b/demos/test_quentin.cpp deleted file mode 100644 index ca05041c..00000000 --- a/demos/test_quentin.cpp +++ /dev/null @@ -1,267 +0,0 @@ -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include - -// int nofriction_two_spheres(int argc, char** argv) -// { -// CLI::App app("Two spheres with no friction"); - -// constexpr std::size_t dim = 2; - -// scopi::scopi_container particles; - -// scopi::sphere s1( -// { -// {1.5, 1.6} -// }, -// 0.5); - -// scopi::sphere s2( -// { -// {3.5, 1.4} -// }, -// 0.5); - -// particles.push_back(s1, -// scopi::property() -// .mass(1) -// .velocity({ -// {1, 0} -// }) -// .moment_inertia(1)); -// particles.push_back(s2, -// scopi::property() -// .mass(1) -// .velocity({ -// {-1, 0} -// }) -// .moment_inertia(1)); - -// double Tf = 14; -// double dt = 0.05; - -// using problem_t = scopi::NoFriction; -// using optim_solver = scopi::OptimGradient; -// using contact_t = scopi::contact_kdtree; -// using vap_t = scopi::vap_fpd; -// scopi::ScopiSolver solver(particles, dt); -// solver.init_options(app); -// CLI11_PARSE(app, argc, argv); -// solver.run(Tf / dt); -// return 0; -// } - -// int nofriction_plan_sphere(int argc, char** argv) -// { -// CLI::App app("Plan/sphere with no friction"); - -// constexpr std::size_t dim = 2; - -// scopi::scopi_container particles; - -// double PI = xt::numeric_constants::PI; -// scopi::plan plan( -// { -// {0., 0.} -// }, -// PI / 2 - PI / 4); -// scopi::sphere sphere( -// { -// {0, 1.5} -// }, -// 0.5); -// particles.push_back(plan, scopi::property().deactivate()); -// particles.push_back(sphere, -// scopi::property().mass(1).moment_inertia(1).velocity({ -// {0, -1} -// })); - -// double Tf = 2.5; -// double dt = 0.05; - -// using problem_t = scopi::NoFriction; -// using optim_solver = scopi::OptimGradient; -// using contact_t = scopi::contact_kdtree; -// using vap_t = scopi::vap_fpd; -// scopi::ScopiSolver solver(particles, dt); -// solver.init_options(app); -// CLI11_PARSE(app, argc, argv); -// solver.run(Tf / dt); -// return 0; -// } - -// int viscous_two_spheres(int argc, char** argv) -// { -// CLI::App app("Two viscous spheres"); - -// constexpr std::size_t dim = 2; - -// scopi::scopi_container particles; - -// scopi::sphere s1( -// { -// {1.5, 1.6} -// }, -// 0.5); - -// scopi::sphere s2( -// { -// {3.5, 1.4} -// }, -// 0.5); - -// particles.push_back(s1, -// scopi::property() -// .mass(1) -// .force({ -// {1, 0} -// }) -// .moment_inertia(1)); -// particles.push_back(s2, -// scopi::property() -// .mass(1) -// .force({ -// {-1, 0} -// }) -// .moment_inertia(1)); - -// double Tf = 5; -// double dt = 0.05; - -// using problem_t = scopi::Viscous; -// using optim_solver = scopi::OptimGradient; -// using contact_t = scopi::contact_kdtree; -// using vap_t = scopi::vap_fpd; -// scopi::ScopiSolver solver(particles, dt); -// solver.init_options(app); -// CLI11_PARSE(app, argc, argv); -// solver.run(Tf / dt); -// return 0; -// } - -int viscous_plan_sphere(int argc, char** argv) -{ - CLI::App app("viscous plan/sphere"); - - constexpr std::size_t dim = 2; - - scopi::scopi_container particles; - - const double PI = xt::numeric_constants::PI; - scopi::plan plan( - { - {0., 0.} - }, - PI / 2); - scopi::sphere sphere( - { - {0, 1.} - }, - 0.5); - particles.push_back(plan, scopi::property().deactivate()); - particles.push_back(sphere, - scopi::property().mass(1).moment_inertia(1).force({ - {0, -1} - })); - - double Tf = 2.5; - double dt = 0.05; - - using problem_t = scopi::Viscous; - using optim_solver = scopi::OptimGradient; - using contact_t = scopi::contact_kdtree; - using vap_t = scopi::vap_fpd; - scopi::ScopiSolver solver(particles, dt); - solver.init_options(app); - CLI11_PARSE(app, argc, argv); - solver.run(Tf / dt); - - Tf = 5.; - particles.f()(1) = {0, 1}; - solver.run(Tf / dt); - - return 0; -} - -int friction_plan_sphere(int argc, char** argv) -{ - CLI::App app("friction plan/sphere"); - - constexpr std::size_t dim = 2; - - scopi::scopi_container particles; - - double PI = xt::numeric_constants::PI; - scopi::plan plan( - { - {0., 0.} - }, - PI / 2 - PI / 6); - - double dd = 1; - double rr = 2; - scopi::sphere sphere( - { - {0, (rr + dd) / std::cos(PI / 6)} - }, - rr); - particles.push_back(plan, scopi::property().deactivate()); - particles.push_back(sphere, - scopi::property() - .mass(1) - .moment_inertia(0.5 * rr * rr) - .force({ - {0, -1} - })); - - double Tf = 10; - double dt = 0.1; - - using problem_t = scopi::Friction; - using optim_solver = scopi::OptimGradient; - using contact_t = scopi::contact_kdtree; - using vap_t = scopi::vap_fpd; - scopi::ScopiSolver solver(particles, dt); - solver.init_options(app); - CLI11_PARSE(app, argc, argv); - solver.run(Tf / dt); - - return 0; -} - -int main(int argc, char** argv) -{ - // static plog::ColorConsoleAppender consoleAppender; - // plog::init(plog::info, &consoleAppender); - - // std::setprecision(15); - // xt::print_options::set_precision(15); - - // nofriction_two_spheres(argc, argv); - // nofriction_plan_sphere(argc, argv); - // viscous_two_spheres(argc, argv); - // viscous_plan_sphere(argc, argv); - friction_plan_sphere(argc, argv); - - return 0; -} diff --git a/demos/three_worms_periodic.cpp b/demos/three_worms_periodic.cpp deleted file mode 100644 index 54c02c0b..00000000 --- a/demos/three_worms_periodic.cpp +++ /dev/null @@ -1,57 +0,0 @@ -#include - -#include -#include -#include - -int main(int argc, char** argv) -{ - plog::init(plog::info, "three_worms.log"); - - CLI::App app("three worms with periodic boundary conditions"); - - constexpr std::size_t dim = 2; - - double dt = .01; - std::size_t total_it = 500; - scopi::scopi_container particles; - auto prop = scopi::property().mass(1.).moment_inertia(0.1); - - double radius = .01; - std::size_t nb_spheres = 6; - std::vector> pos1(nb_spheres), pos2(nb_spheres), pos3(nb_spheres); - std::vector q(nb_spheres); - for (std::size_t i = 0; i < pos1.size(); ++i) - { - pos1[i](0) = 0.42 + 0.52 - i * 2 * radius; - pos1[i][1] = .495; - pos2[i][0] = 0.58 - 0.48 + i * 2 * radius; - pos2[i][1] = .505; - pos3[i][0] = .5 - 0.48; - pos3[i][1] = 0.5 - i * 2 * radius; - - // pos1[i](0) = 0.42 - i*2*radius; - // pos1[i][1] = .495; - // pos2[i][0] = 0.58 + i*2*radius; - // pos2[i][1] = .505; - // pos3[i][0] = .5 ; - // pos3[i][1] = 0.5 - i*2*radius; - q[i] = scopi::quaternion(0); - } - - scopi::worm w1(pos1, q, radius, nb_spheres); - scopi::worm w2(pos2, q, radius, nb_spheres); - scopi::worm w3(pos3, q, radius, nb_spheres); - particles.push_back(w1, prop.desired_velocity({.5, 0.})); - particles.push_back(w2, prop.desired_velocity({-.5, 0.})); - particles.push_back(w3, prop.desired_velocity({0., 0.})); - - auto domain = scopi::BoxDomain({0, 0}, {1, 1}).with_periodicity(0); - scopi::ScopiSolver> solver(domain, particles, dt); - - solver.init_options(app); - CLI11_PARSE(app, argc, argv); - solver.run(total_it); - - return 0; -} diff --git a/demos/two_ellipsoids.cpp b/demos/two_ellipsoids.cpp index e5ccd6fe..e6aae54e 100644 --- a/demos/two_ellipsoids.cpp +++ b/demos/two_ellipsoids.cpp @@ -1,15 +1,12 @@ #include -#include -#include - #include -#include +#include #include int main() { - plog::init(plog::info, "two_ellipsoids.log"); + scopi::initialize("Two ellipsoids"); constexpr std::size_t dim = 2; double PI = xt::numeric_constants::PI; @@ -48,7 +45,7 @@ int main() particles.push_back(s1, prop1); particles.push_back(s2, prop2); - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + solver.run(dt, total_it); return 0; } diff --git a/demos/two_ellipsoids_asymmetrical_2d.cpp b/demos/two_ellipsoids_asymmetrical_2d.cpp deleted file mode 100644 index 8d91df50..00000000 --- a/demos/two_ellipsoids_asymmetrical_2d.cpp +++ /dev/null @@ -1,48 +0,0 @@ -#include -#include -#include - -int main() -{ - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - double dt = .005; - std::size_t total_it = 1000; - scopi::scopi_container particles; - - scopi::superellipsoid s1( - { - {-0.2, 0.} - }, - {scopi::quaternion(PI / 4)}, - {{.1, .05}}, - 1); - scopi::superellipsoid s2( - { - {0.2, 0.02} - }, - {scopi::quaternion(-PI / 4)}, - {{.1, .05}}, - 1); - - auto prop1 = scopi::property() - .desired_velocity({ - {0.25, 0} - }) - .mass(1.) - .moment_inertia(0.1); - auto prop2 = scopi::property() - .desired_velocity({ - {-0.25, 0} - }) - .mass(1.) - .moment_inertia(0.1); - - particles.push_back(s1, prop1); - particles.push_back(s2, prop2); - - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); - - return 0; -} diff --git a/demos/two_ellipsoids_asymmetrical_3d.cpp b/demos/two_ellipsoids_asymmetrical_3d.cpp deleted file mode 100644 index f8e793a7..00000000 --- a/demos/two_ellipsoids_asymmetrical_3d.cpp +++ /dev/null @@ -1,48 +0,0 @@ -#include -#include -#include - -int main() -{ - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - double dt = .005; - std::size_t total_it = 1000; - scopi::scopi_container particles; - - scopi::superellipsoid s1( - { - {-0.2, 0., 0.} - }, - {scopi::quaternion(PI / 4)}, - {{.1, .05, .05}}, - {{1, 1}}); - scopi::superellipsoid s2( - { - {0.2, 0.02, 0.} - }, - {scopi::quaternion(-PI / 4)}, - {{.1, .05, .05}}, - {{1, 1}}); - - auto prop1 = scopi::property() - .desired_velocity({ - {0.25, 0, 0} - }) - .mass(1.) - .moment_inertia(0.1); - auto prop2 = scopi::property() - .desired_velocity({ - {-0.25, 0, 0} - }) - .mass(1.) - .moment_inertia(0.1); - - particles.push_back(s1, prop1); - particles.push_back(s2, prop2); - - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); - - return 0; -} diff --git a/demos/two_ellipsoids_spheres_2d.cpp b/demos/two_ellipsoids_spheres_2d.cpp deleted file mode 100644 index 88ff402d..00000000 --- a/demos/two_ellipsoids_spheres_2d.cpp +++ /dev/null @@ -1,48 +0,0 @@ -#include -#include -#include - -int main() -{ - constexpr std::size_t dim = 2; - double PI = xt::numeric_constants::PI; - double dt = .005; - std::size_t total_it = 1000; - scopi::scopi_container particles; - - scopi::superellipsoid s1( - { - {-0.2, 0.} - }, - {scopi::quaternion(PI / 4)}, - {{.1, .1}}, - 1); - scopi::superellipsoid s2( - { - {0.2, 0.} - }, - {scopi::quaternion(-PI / 4)}, - {{.1, .1}}, - 1); - - auto prop1 = scopi::property() - .desired_velocity({ - {0.25, 0} - }) - .mass(1.) - .moment_inertia(0.1); - auto prop2 = scopi::property() - .desired_velocity({ - {-0.25, 0} - }) - .mass(1.) - .moment_inertia(0.1); - - particles.push_back(s1, prop1); - particles.push_back(s2, prop2); - - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); - - return 0; -} diff --git a/demos/two_ellipsoids_spheres_3d.cpp b/demos/two_ellipsoids_spheres_3d.cpp deleted file mode 100644 index dcacb963..00000000 --- a/demos/two_ellipsoids_spheres_3d.cpp +++ /dev/null @@ -1,48 +0,0 @@ -#include -#include -#include - -int main() -{ - constexpr std::size_t dim = 3; - double PI = xt::numeric_constants::PI; - double dt = .005; - std::size_t total_it = 1000; - scopi::scopi_container particles; - - scopi::superellipsoid s1( - { - {-0.2, 0., 0.} - }, - {scopi::quaternion(PI / 4)}, - {{.1, .1, .1}}, - {{1, 1}}); - scopi::superellipsoid s2( - { - {0.2, 0., 0.} - }, - {scopi::quaternion(-PI / 4)}, - {{.1, .1, .1}}, - {{1, 1}}); - - auto prop1 = scopi::property() - .desired_velocity({ - {0.25, 0, 0} - }) - .mass(1.) - .moment_inertia(0.1); - auto prop2 = scopi::property() - .desired_velocity({ - {-0.25, 0, 0} - }) - .mass(1.) - .moment_inertia(0.1); - - particles.push_back(s1, prop1); - particles.push_back(s2, prop2); - - scopi::ScopiSolver solver(particles, dt); - solver.run(total_it); - - return 0; -} diff --git a/demos/two_spheres.cpp b/demos/two_spheres.cpp index 17e9bf33..0e30bac2 100644 --- a/demos/two_spheres.cpp +++ b/demos/two_spheres.cpp @@ -1,20 +1,13 @@ +#include #include -#include +#include #include -#include - -// #include int main(int argc, char** argv) { - plog::init(plog::info, "two_spheres.log"); - CLI::App app("two spheres"); + scopi::initialize("Two spheres simulation"); constexpr std::size_t dim = 2; - double dt = .005; - std::size_t total_it = 1000; - scopi::scopi_container particles; - scopi::sphere s1( { {-0.2, -0.05} @@ -25,6 +18,11 @@ int main(int argc, char** argv) {0.2, 0.05} }, 0.1); + + scopi::scopi_container particles; + scopi::ScopiSolver solver(particles); + SCOPI_PARSE(argc, argv); + particles.push_back(s1, scopi::property() .desired_velocity({ @@ -40,10 +38,10 @@ int main(int argc, char** argv) .mass(1.) .moment_inertia(0.1)); - scopi::ScopiSolver solver(particles, dt); - solver.init_options(app); - CLI11_PARSE(app, argc, argv); - solver.run(total_it); + double dt = 0.005; + + std::size_t total_it = 100; + solver.run(dt, total_it); return 0; } diff --git a/demos/two_spheres_periodic.cpp b/demos/two_spheres_periodic.cpp deleted file mode 100644 index de396165..00000000 --- a/demos/two_spheres_periodic.cpp +++ /dev/null @@ -1,58 +0,0 @@ -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -int main(int argc, char** argv) -{ - plog::init(plog::info, "two_spheres_periodic.log", 0, 0); - - CLI::App app("two spheres with periodic boundary conditions"); - constexpr std::size_t dim = 2; - double dt = .005; - std::size_t total_it = 1000; - - double vel = -0.25; - - scopi::scopi_container particles; - - scopi::sphere s1( - { - {0.4, 0.45} - }, - 0.1); - scopi::sphere s2( - { - {0.6, 0.55} - }, - 0.1); - particles.push_back(s1, - scopi::property() - .desired_velocity({ - {vel, 0} - }) - .mass(1.) - .moment_inertia(0.1)); - particles.push_back(s2, - scopi::property() - .desired_velocity({ - {-vel, 0} - }) - .mass(1.) - .moment_inertia(0.1)); - - auto domain = scopi::BoxDomain({0, 0}, {1, 1}).with_periodicity(0); - scopi::ScopiSolver solver(domain, particles, dt); - solver.init_options(app); - CLI11_PARSE(app, argc, argv); - solver.run(total_it); - - return 0; -} diff --git a/demos/two_spheres_viscosity.cpp b/demos/two_spheres_viscosity.cpp deleted file mode 100644 index a3d48eb2..00000000 --- a/demos/two_spheres_viscosity.cpp +++ /dev/null @@ -1,75 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include - -int main() -{ - plog::init(plog::info, "two_spheres_viscosity.log"); - - constexpr std::size_t dim = 2; - - const double dt = .001; - const std::size_t total_it = 3000; - - const double r = 0.1; - scopi::scopi_container particles; - auto prop = scopi::property().mass(1.).moment_inertia(1. * r * r / 2.); - - scopi::sphere s1( - { - {0., 0.} - }, - r); - scopi::sphere s2( - { - {1., 0.05} - }, - r); - scopi::sphere s3( - { - {0., 0.3} - }, - r); - scopi::sphere s4( - { - {1., 0.35} - }, - r); - // particles.push_back(s1, prop.velocity({1., 0.})); - // particles.push_back(s1, scopi::property().deactivate()); - // particles.push_back(s2, prop.velocity({-1., 0.})); - particles.push_back(s1, - prop.force({ - {1, 0} - })); - particles.push_back(s2, - prop.force({ - {-1, 0} - })); - particles.push_back(s3, - prop.force({ - {1, 0} - })); - particles.push_back(s4, - prop.force({ - {-1, 0} - })); - - scopi::ScopiSolver>, scopi::contact_kdtree, scopi::vap_fpd> solver(particles, dt); - auto params = solver.get_params(); - params.optim_params.change_default_tol_mosek = false; - params.problem_params.mu = 0.1; - params.contact_params.dmax = 0.11; - params.problem_params.gamma_min = -0.7; - // params.problem_params.tol = 10^(-8); - params.solver_params.output_frequency = 1; - params.solver_params.path = "_two_spheres"; - params.solver_params.filename = "two_spheres"; - solver.run(total_it); - - return 0; -} diff --git a/demos/two_worms.cpp b/demos/two_worms.cpp index 420ec901..a777c6c5 100644 --- a/demos/two_worms.cpp +++ b/demos/two_worms.cpp @@ -1,18 +1,11 @@ -#include -#include - -#include +#include #include -#include +#include #include -#include -#include int main(int argc, char** argv) { - plog::init(plog::error, "two_worms.log"); - - CLI::App app("two spheres with periodic boundary conditions"); + scopi::initialize("Two worms simulation"); constexpr std::size_t dim = 2; double dt = .005; @@ -57,10 +50,10 @@ int main(int argc, char** argv) particles.push_back(w1, prop.desired_velocity({-1., 0.})); particles.push_back(w2, prop.desired_velocity({1., 0.})); - scopi::ScopiSolver, scopi::contact_brute_force> solver(particles, dt); - solver.init_options(app); - CLI11_PARSE(app, argc, argv); - solver.run(total_it); + scopi::ScopiSolver solver(particles); + SCOPI_PARSE(argc, argv); + + solver.run(dt, total_it); return 0; } diff --git a/demos/viscosity.cpp b/demos/viscosity.cpp deleted file mode 100644 index 81f2de11..00000000 --- a/demos/viscosity.cpp +++ /dev/null @@ -1,56 +0,0 @@ -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -int main() -{ - plog::init(plog::info, "viscosity.log"); - - constexpr std::size_t dim = 2; - using solver_t = scopi::ScopiSolver>, scopi::contact_kdtree, scopi::vap_fpd>; - double PI = xt::numeric_constants::PI; - - double radius = 1.; - double g = 1.; - double h = 1.5 * radius; - double alpha = PI / 4.; - auto prop = scopi::property().mass(1.).moment_inertia(1. * radius * radius / 2.); - - double dt = 0.05; - std::size_t total_it = 60; - scopi::scopi_container particles; - scopi::plan p( - { - {0., 0.} - }, - PI / 2.); - scopi::sphere s( - { - {0., 2 * h} - }, - radius); - particles.push_back(p, scopi::property().deactivate()); - particles.push_back(s, - prop.force({ - {g * std::cos(alpha), -g * std::sin(alpha)} - })); - solver_t solver(particles, dt); - auto params = solver.get_params(); - params.contact_params.dmax = 2; - solver.run(total_it); - particles.f()(1)(1) *= -1.; - solver.run(3 * total_it, total_it); - particles.f()(1)(1) *= -1.; - solver.run(8 * total_it, 3 * total_it); - - return 0; -} diff --git a/doc/source/api/objects/types/plan.rst b/doc/source/api/objects/types/plan.rst index 7f4b844b..978423b4 100644 --- a/doc/source/api/objects/types/plan.rst +++ b/doc/source/api/objects/types/plan.rst @@ -1,6 +1,6 @@ -plan class +plane class ========== -.. doxygenclass:: scopi::plan +.. doxygenclass:: scopi::plane :project: scopi :members: diff --git a/doc/source/api/problems/DryWithFriction.rst b/doc/source/api/problems/DryWithFriction.rst index 917f6080..3ae36dfb 100644 --- a/doc/source/api/problems/DryWithFriction.rst +++ b/doc/source/api/problems/DryWithFriction.rst @@ -12,7 +12,7 @@ analytical_solution_sphere_plan functions .. doxygenfunction:: scopi::analytical_solution_sphere_plan :project: scopi -.. doxygenfunction:: scopi::analytical_solution_sphere_plan_velocity +.. doxygenfunction:: scopi::analytical_solution_sphere_plane_velocity :project: scopi ProblemParams class diff --git a/doc/source/index.rst b/doc/source/index.rst index 71233724..147f67c3 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -59,7 +59,7 @@ Objects api/objects/neighbor api/objects/types/sphere api/objects/types/superellipsoid - api/objects/types/plan + api/objects/types/plane api/objects/types/worm Problems diff --git a/include/scopi/contact/contact_kdtree.hpp b/include/scopi/contact/contact_kdtree.hpp index 9df0d042..02bd6978 100644 --- a/include/scopi/contact/contact_kdtree.hpp +++ b/include/scopi/contact/contact_kdtree.hpp @@ -54,7 +54,7 @@ namespace scopi * @brief Kd-tree radius. * * For two particles \c i and \c j, compute the exact distance only if the squared distance between a point in \c i and a point in - * \c j is less than \c kdtree_radius. For a sphere or a superellipsoid, this point is the center. For a plan, it is the point used + * \c j is less than \c kdtree_radius. For a sphere or a superellipsoid, this point is the center. For a plane, it is the point used * to construct it. \note \c kd_tree_radius > 0 */ double kd_tree_radius{17.}; diff --git a/include/scopi/container.hpp b/include/scopi/container.hpp index 981fc4fb..760af864 100644 --- a/include/scopi/container.hpp +++ b/include/scopi/container.hpp @@ -32,7 +32,7 @@ namespace scopi * of the container. * * In the following, "particle" means a base object (sphere, superellipsoid - * or plan) and an "object" can be a more complex object, such as a worm. + * or plane) and an "object" can be a more complex object, such as a worm. * Particles are objects. * * @tparam dim Dimension (2 or 3). @@ -75,7 +75,7 @@ namespace scopi /** * @brief Reconstructs an object. * - * An object can be a sphere, a superellipsoid, a plan, a worm,... + * An object can be a sphere, a superellipsoid, a plane, a worm,... * * @param i Index of the object. * diff --git a/include/scopi/objects/methods/closest_points.hpp b/include/scopi/objects/methods/closest_points.hpp index fdcf1afe..16fb069f 100644 --- a/include/scopi/objects/methods/closest_points.hpp +++ b/include/scopi/objects/methods/closest_points.hpp @@ -13,7 +13,7 @@ #include "../dispatch.hpp" #include "../neighbor.hpp" -#include "../types/plan.hpp" +#include "../types/plane.hpp" #include "../types/segment.hpp" #include "../types/sphere.hpp" #include "../types/superellipsoid.hpp" @@ -732,7 +732,7 @@ namespace scopi return neigh; } - // PLAN - PLAN + // PLANE - PLANE /** * @brief Neighbor between two planes. * @@ -747,12 +747,12 @@ namespace scopi * @return Neighbor struct for contact between plane \c i and plane \c j. */ template - auto closest_points(const plan& /* p1 */, const plan& /* p2 */) + auto closest_points(const plane& /* p1 */, const plane& /* p2 */) { return neighbor(); } - // SPHERE - PLAN + // SPHERE - PLANE /** * @brief Neighbor between a sphere and a plane. * @@ -766,27 +766,27 @@ namespace scopi * @return Neighbor struct for contact between sphere \c i and plane \c j. */ template - auto closest_points(const sphere& s, const plan& p) + auto closest_points(const sphere& s, const plane& p) { - // std::cout << "closest_points : SPHERE - PLAN" << std::endl; + // std::cout << "closest_points : SPHERE - PLANE" << std::endl; auto s_pos = s.pos(0); auto p_pos = p.pos(0); auto normal = p.normal(); // plan2sphs.n - auto plan_to_sphere = xt::eval(xt::linalg::dot(s_pos - p_pos, normal)); - auto xtsign = xt::sign(plan_to_sphere); + auto plane_to_sphere = xt::eval(xt::linalg::dot(s_pos - p_pos, normal)); + auto xtsign = xt::sign(plane_to_sphere); neighbor neigh; neigh.pi = s_pos - xtsign * s.radius() * normal; - neigh.pj = s_pos - plan_to_sphere * normal; + neigh.pj = s_pos - plane_to_sphere * normal; neigh.nij = xtsign * normal; neigh.dij = xt::linalg::dot(neigh.pi - neigh.pj, neigh.nij)[0]; return neigh; } - // PLAN - SPHERE + // PLANE - SPHERE /** * @brief Neighbor between a plane and a sphere. * @@ -800,7 +800,7 @@ namespace scopi * @return Neighbor struct for contact between plane \c i and sphere \c j. */ template - auto closest_points(const plan& p, const sphere& s) + auto closest_points(const plane& p, const sphere& s) { auto neigh = closest_points(s, p); neigh.nij *= -1.; @@ -902,16 +902,16 @@ namespace scopi } template - auto closest_points(const plan&, const segment&) + auto closest_points(const plane&, const segment&) { - std::cerr << "closest_points for plan/segment must be implemented" << std::endl; + std::cerr << "closest_points for plane/segment must be implemented" << std::endl; return neighbor(); } template - auto closest_points(const segment&, const plan&) + auto closest_points(const segment&, const plane&) { - std::cerr << "closest_points for segment/plan must be implemented" << std::endl; + std::cerr << "closest_points for segment/plane must be implemented" << std::endl; return neighbor(); } @@ -1463,7 +1463,7 @@ namespace scopi return neigh; } - // SUPERELLIPSOID 3D - PLAN 3D + // SUPERELLIPSOID 3D - PLANE 3D /** * @brief Neighbor between a superellipsoid and a plane in 3D. * @@ -1476,9 +1476,9 @@ namespace scopi * @return Neighbor struct for contact between superellipsod \c i and plane \c j. */ template - auto closest_points(const superellipsoid<3, owner> s1, const plan<3, owner> p2) + auto closest_points(const superellipsoid<3, owner> s1, const plane<3, owner> p2) { - std::cout << "closest_points : SUPERELLIPSOID 3D - PLAN 3D" << std::endl; + std::cout << "closest_points : SUPERELLIPSOID 3D - PLANE 3D" << std::endl; double pi = 4 * std::atan(1); neighbor<3, problem_t> neigh; // Pour déterminer de quel cote on est @@ -1745,7 +1745,7 @@ namespace scopi return neigh; } - // PLAN 3D - SUPERELLIPSOID 3D + // PLANE 3D - SUPERELLIPSOID 3D /** * @brief Neighbor between a plane and a superellipsoid in 3D. * @@ -1759,7 +1759,7 @@ namespace scopi * @return Neighbor struct for contact between plane \c i and superellipsoid \c j. */ template - auto closest_points(const plan<3, owner> p2, const superellipsoid<3, owner> s1) + auto closest_points(const plane<3, owner> p2, const superellipsoid<3, owner> s1) { auto neigh = closest_points(s1, p2); neigh.nij *= -1.; @@ -1780,7 +1780,7 @@ namespace scopi */ // SUPERELLIPSOID 2D - DROITE 2D template - auto closest_points(const superellipsoid<2, owner> s1, const plan<2, owner> d2) + auto closest_points(const superellipsoid<2, owner> s1, const plane<2, owner> d2) { // std::cout << "closest_points : SUPERELLIPSOID 2D - DROITE 2D" << std::endl; double pi = 4 * std::atan(1); @@ -1949,7 +1949,7 @@ namespace scopi */ // DROITE 2D - SUPERELLIPSOID 2D template - auto closest_points(const plan<2, owner> d2, const superellipsoid<2, owner> s1) + auto closest_points(const plane<2, owner> d2, const superellipsoid<2, owner> s1) { auto neigh = closest_points(s1, d2); neigh.nij *= -1.; @@ -2038,7 +2038,7 @@ namespace scopi using closest_points_dispatcher = double_static_dispatcher< closest_points_functor, const object, - mpl::vector, const superellipsoid, const plan, const segment>, + mpl::vector, const superellipsoid, const plane, const segment>, typename closest_points_functor::return_type, antisymmetric_dispatch>; } diff --git a/include/scopi/objects/methods/distance.hpp b/include/scopi/objects/methods/distance.hpp index 390c6302..440cf3d1 100644 --- a/include/scopi/objects/methods/distance.hpp +++ b/include/scopi/objects/methods/distance.hpp @@ -4,7 +4,7 @@ #include "../dispatch.hpp" #include "../types/globule.hpp" -#include "../types/plan.hpp" +#include "../types/plane.hpp" #include "../types/sphere.hpp" #include "../types/superellipsoid.hpp" @@ -26,11 +26,11 @@ namespace scopi return xt::linalg::norm(s1.pos() - s2.pos()) - s1.radius() - s2.radius(); } - // PLAN - PLAN + // PLANE - PLANE template - double distance(const plan, const plan) + double distance(const plane, const plane) { - std::cout << "distance : PLAN - PLAN" << std::endl; + std::cout << "distance : PLANE - PLANE" << std::endl; return 20; } @@ -58,11 +58,11 @@ namespace scopi return 50; } - // SPHERE - PLAN + // SPHERE - PLANE template - double distance(const sphere, const plan) + double distance(const sphere, const plane) { - std::cout << "distance : SPHERE - PLAN" << std::endl; + std::cout << "distance : SPHERE - PLANE" << std::endl; return 60; } @@ -74,19 +74,19 @@ namespace scopi return 70; } - // SUPERELLIPSOID - PLAN + // SUPERELLIPSOID - PLANE template - double distance(const superellipsoid, const plan) + double distance(const superellipsoid, const plane) { - std::cout << "distance : SUPERELLIPSOID - PLAN" << std::endl; + std::cout << "distance : SUPERELLIPSOID - PLANE" << std::endl; return 80; } - // GLOBULE - PLAN + // GLOBULE - PLANE template - double distance(const globule, const plan) + double distance(const globule, const plane) { - std::cout << "distance : GLOBULE - PLAN" << std::endl; + std::cout << "distance : GLOBULE - PLANE" << std::endl; return 90; } @@ -111,7 +111,7 @@ namespace scopi using distance_dispatcher = double_static_dispatcher< distance_functor, const object, - mpl::vector, const superellipsoid, const globule, const plan>, + mpl::vector, const superellipsoid, const globule, const plane>, typename distance_functor::return_type, symmetric_dispatch>; } diff --git a/include/scopi/objects/methods/select.hpp b/include/scopi/objects/methods/select.hpp index 12f2ff1e..7e86da05 100644 --- a/include/scopi/objects/methods/select.hpp +++ b/include/scopi/objects/methods/select.hpp @@ -13,7 +13,7 @@ #include "../dispatch.hpp" #include "../neighbor.hpp" -#include "../types/plan.hpp" +#include "../types/plane.hpp" #include "../types/segment.hpp" #include "../types/sphere.hpp" #include "../types/superellipsoid.hpp" @@ -50,11 +50,11 @@ namespace scopi return std::make_unique>(s); } - // PLAN + // PLANE template - std::unique_ptr> select_object(const plan& s, const std::size_t) + std::unique_ptr> select_object(const plane& s, const std::size_t) { - return std::make_unique>(s); + return std::make_unique>(s); } // SEGMENT @@ -92,7 +92,7 @@ namespace scopi using select_object_dispatcher = double_static_dispatcher< select_object_functor, const object, - mpl::vector, const superellipsoid, const worm, const plan, const segment>, + mpl::vector, const superellipsoid, const worm, const plane, const segment>, typename select_object_functor::return_type, antisymmetric_dispatch, const index, diff --git a/include/scopi/objects/methods/write_objects.hpp b/include/scopi/objects/methods/write_objects.hpp index e162c32c..bdf1efe8 100644 --- a/include/scopi/objects/methods/write_objects.hpp +++ b/include/scopi/objects/methods/write_objects.hpp @@ -12,7 +12,7 @@ #include "../dispatch.hpp" #include "../neighbor.hpp" -#include "../types/plan.hpp" +#include "../types/plane.hpp" #include "../types/segment.hpp" #include "../types/sphere.hpp" #include "../types/superellipsoid.hpp" @@ -104,7 +104,7 @@ namespace scopi // return ssss; } - // PLAN + // PLANE /** * @brief * @brief Write the elements of a plane in json format. @@ -115,11 +115,11 @@ namespace scopi * @return nlohmann json object. */ template - nl::json write_objects(const plan& p, [[maybe_unused]] std::size_t id) + nl::json write_objects(const plane& p, [[maybe_unused]] std::size_t id) { nl::json object; - object["type"] = "plan"; + object["type"] = "plane"; object["id"] = id; object["position"] = p.pos(); object["normal"] = p.normal(); @@ -128,7 +128,7 @@ namespace scopi return object; - // std::cout << "write_objects : PLAN" << std::endl; + // std::cout << "write_objects : PLANE" << std::endl; // std::stringstream ss; // std::regex xp1("\\{|\\}"); // std::regex xp2("\\.,"); @@ -257,6 +257,6 @@ namespace scopi using write_objects_dispatcher = unit_static_dispatcher< write_objects_functor, const object, - mpl::vector, const superellipsoid, const worm, const plan, const segment>, + mpl::vector, const superellipsoid, const worm, const plane, const segment>, typename write_objects_functor::return_type>; } diff --git a/include/scopi/objects/types/plan.hpp b/include/scopi/objects/types/plane.hpp similarity index 81% rename from include/scopi/objects/types/plan.hpp rename to include/scopi/objects/types/plane.hpp index 7296883e..df558e99 100644 --- a/include/scopi/objects/types/plan.hpp +++ b/include/scopi/objects/types/plane.hpp @@ -13,17 +13,17 @@ namespace scopi { ///////////////////// - // plan definition // + // plane definition // ///////////////////// /** - * @class plan + * @class plane * @brief Plane. * * @tparam dim Dimension (2 or 3). * @tparam owner */ template - class plan : public object + class plane : public object { public: @@ -46,14 +46,14 @@ namespace scopi * @param pos [in] Position of a point in the plane. * @param angle [in] Angle of the rotation. */ - explicit plan(position_type pos, double angle = 0); + explicit plane(position_type pos, double angle = 0); /** * @brief Constructor with given rotation. * * @param pos [in] Position of a point in the plane. * @param q [in] Quaternion describing the rotation of the plane. */ - explicit plan(position_type pos, quaternion_type q); + explicit plane(position_type pos, quaternion_type q); /** * @brief Get the rotation matrix of the plane. @@ -117,63 +117,63 @@ namespace scopi }; ///////////////////////// - // plan implementation // + // plane implementation // ///////////////////////// template - plan::plan(position_type pos, double angle) + plane::plane(position_type pos, double angle) : base_type(pos, {quaternion(angle)}, 1) { create_hash(); } template - auto plan::rotation() const + auto plane::rotation() const { return rotation_matrix(this->q()); } template - plan::plan(position_type pos, quaternion_type q) + plane::plane(position_type pos, quaternion_type q) : base_type(pos, q, 1) { create_hash(); } template - auto plan::normal() const + auto plane::normal() const { auto matrix = rotation_matrix(this->q(0)); return xt::eval(xt::view(matrix, xt::all(), 0)); } template - std::unique_ptr> plan::construct() const + std::unique_ptr> plane::construct() const { - return make_object_constructor>(); + return make_object_constructor>(); } template - void plan::print() const + void plane::print() const { - std::cout << "plan<" << dim << ">()\n"; + std::cout << "plane<" << dim << ">()\n"; } template - std::size_t plan::hash() const + std::size_t plane::hash() const { return m_hash; } template - void plan::create_hash() + void plane::create_hash() { std::stringstream ss; - ss << "plan<" << dim << ">()"; + ss << "plane<" << dim << ">()"; m_hash = std::hash{}(ss.str()); } template - auto plan::point(double b) const + auto plane::point(double b) const { xt::xtensor_fixed> pt; pt(0) = 0; @@ -182,7 +182,7 @@ namespace scopi } template - auto plan::point(double a, double b) const + auto plane::point(double a, double b) const { xt::xtensor_fixed> pt; pt(0) = 0; diff --git a/include/scopi/solver.hpp b/include/scopi/solver.hpp index bf08009d..e55fc443 100644 --- a/include/scopi/solver.hpp +++ b/include/scopi/solver.hpp @@ -114,7 +114,7 @@ namespace scopi * @param dt Time step. It is fixed during the simulation. * @param params Parameters for the different steps of the algorithm. */ - ScopiSolver(const BoxDomain& box, scopi_container& particles, double dt); + explicit ScopiSolver(const BoxDomain& box, scopi_container& particles); /** * @brief Constructor. @@ -123,7 +123,7 @@ namespace scopi * @param dt Time step. It is fixed during the simulation. * @param params Parameters for the different steps of the algorithm. */ - ScopiSolver(scopi_container& particles, double dt); + explicit ScopiSolver(scopi_container& particles); void init_options(); @@ -133,7 +133,7 @@ namespace scopi * @param total_it [in] Total number of iterations to perform. * @param initial_iter [in] Initial index of iteration. Used for restart or to change external parameters. */ - void run(std::size_t total_it, std::size_t initial_iter = 0); + void run(double dt, std::size_t total_it, std::size_t initial_iter = 0); /** * @brief Return the current contacts of the simulation. @@ -149,6 +149,8 @@ namespace scopi private: + void set_timestep(double dt); + /** * @brief Move obstacles (particles with an imposed velocity). */ @@ -194,7 +196,7 @@ namespace scopi /** * @brief Time step, fixed during the simulation. */ - double m_dt; + double m_dt = 0; optim_solver_t m_optim_solver; contact_method_t m_contact_method; @@ -205,25 +207,22 @@ namespace scopi template class contact_method_t, class vap_t> ScopiSolver::ScopiSolver(const BoxDomain& box, - scopi_container& particles, - double dt) + scopi_container& particles) : m_box(box) , m_particles(particles) - , m_dt(dt) - , m_optim_solver(particles.nb_active(), dt, particles) + , m_optim_solver() , m_contact_method() - , m_vap(particles.nb_active(), particles.nb_inactive(), particles.size(), dt) + , m_vap() { init_options(); } template class contact_method_t, class vap_t> - ScopiSolver::ScopiSolver(scopi_container& particles, double dt) + ScopiSolver::ScopiSolver(scopi_container& particles) : m_particles(particles) - , m_dt(dt) - , m_optim_solver(particles.nb_active(), dt, particles) + , m_optim_solver() , m_contact_method() - , m_vap(particles.nb_active(), particles.nb_inactive(), particles.size(), dt) + , m_vap() { init_options(); } @@ -259,10 +258,19 @@ namespace scopi } template class contact_method_t, class vap_t> - void ScopiSolver::run(std::size_t total_it, std::size_t initial_iter) + void ScopiSolver::set_timestep(double dt) + { + m_dt = dt; + } + + template class contact_method_t, class vap_t> + void ScopiSolver::run(double dt, std::size_t total_it, std::size_t initial_iter) { // Time Loop write_output_files(m_old_contacts, initial_iter); + set_timestep(dt); + m_optim_solver.set_timestep(m_dt); + for (std::size_t nite = initial_iter; nite < total_it; ++nite) { PLOG_INFO << "\n\n------------------- Time iteration ----------------> " << nite; @@ -275,7 +283,7 @@ namespace scopi transfer(m_old_contacts, contacts); } - m_vap.set_a_priori_velocity(m_particles, contacts); + m_vap.set_a_priori_velocity(m_dt, m_particles, contacts); m_optim_solver.extra_steps_before_solve(contacts); while (m_optim_solver.should_solve()) { @@ -422,7 +430,16 @@ namespace scopi #pragma omp parallel for for (std::size_t i = 0; i < m_particles.nb_active(); ++i) { - xt::xtensor_fixed> w({0, 0, m_particles.omega()(i + active_offset)}); + xt::xtensor_fixed> w; + if constexpr (dim == 2) + { + w = {0, 0, m_particles.omega()(i + active_offset)}; + } + else + { + w = m_particles.omega()(i + active_offset); + } + double normw = xt::linalg::norm(w); if (normw == 0) { diff --git a/include/scopi/solvers/OptimGradient.hpp b/include/scopi/solvers/OptimGradient.hpp index 0907b72d..9b65bd8b 100644 --- a/include/scopi/solvers/OptimGradient.hpp +++ b/include/scopi/solvers/OptimGradient.hpp @@ -61,11 +61,9 @@ namespace scopi using method_t = Method; using params_t = typename method_t::params_t; - template - OptimGradient(std::size_t, double dt, const scopi_container&) - : m_dt(dt) - , m_method(method_t()) + void set_timestep(double dt) { + m_dt = dt; } void init_options() diff --git a/include/scopi/vap/base.hpp b/include/scopi/vap/base.hpp index 7c86d7c8..01bd38a4 100644 --- a/include/scopi/vap/base.hpp +++ b/include/scopi/vap/base.hpp @@ -28,53 +28,21 @@ namespace scopi * @param contacts [in] Array of contacts. */ template - void set_a_priori_velocity(scopi_container& particles, const Contacts& contacts); - /** - * @brief Constructor. - * - * @param Nactive Number of active particles. - * @param active_ptr Index of the first active particle. - * @param dt Time step. - */ - vap_base(std::size_t Nactive, std::size_t active_ptr, double dt); + void set_a_priori_velocity(double dt, scopi_container& particles, const Contacts& contacts); params_t& get_params(); - protected: - - /** - * @brief Number of active particles. - */ - std::size_t m_Nactive; - /** - * @brief Index of the first active particle. - */ - std::size_t m_active_ptr; - /** - * @brief Time step. - */ - double m_dt; - private: params_t m_params; }; - template - vap_base::vap_base(std::size_t Nactive, std::size_t active_ptr, double dt) - : m_Nactive(Nactive) - , m_active_ptr(active_ptr) - , m_dt(dt) - , m_params() - { - } - template template - void vap_base::set_a_priori_velocity(scopi_container& particles, const Contacts& contacts) + void vap_base::set_a_priori_velocity(double dt, scopi_container& particles, const Contacts& contacts) { tic(); - this->derived_cast().set_a_priori_velocity_impl(particles, contacts); + this->derived_cast().set_a_priori_velocity_impl(dt, particles, contacts); auto duration = toc(); PLOG_INFO << "----> CPUTIME : set vap = " << duration; } diff --git a/include/scopi/vap/vap_fixed.hpp b/include/scopi/vap/vap_fixed.hpp index 07ff6998..fce12ecb 100644 --- a/include/scopi/vap/vap_fixed.hpp +++ b/include/scopi/vap/vap_fixed.hpp @@ -39,24 +39,15 @@ namespace scopi * @param contacts [in] Array of contacts. */ template - void set_a_priori_velocity_impl(scopi_container& particles, const Contacts& contacts); - - /** - * @brief Constructor. - * - * @param Nactive Number of active particles. - * @param active_ptr Index of the first active particle. - * @param nb_parts Number of objects. - * @param dt Time step. - * @param params Parameters (for compatibility). - */ - vap_fixed(std::size_t Nactive, std::size_t active_ptr, std::size_t nb_parts, double dt); + void set_a_priori_velocity_impl(double dt, scopi_container& particles, const Contacts& contacts); }; template - void vap_fixed::set_a_priori_velocity_impl(scopi_container& particles, const Contacts&) + void vap_fixed::set_a_priori_velocity_impl(double, scopi_container& particles, const Contacts&) { - for (std::size_t i = m_active_ptr; i < m_active_ptr + m_Nactive; ++i) + auto active_ptr = particles.nb_inactive(); + auto nb_active = particles.nb_active(); + for (std::size_t i = active_ptr; i < active_ptr + nb_active; ++i) { particles.v()(i) = particles.vd()(i); particles.omega()(i) = particles.desired_omega()(i); diff --git a/include/scopi/vap/vap_fpd.hpp b/include/scopi/vap/vap_fpd.hpp index 76b6d6f5..1967e81e 100644 --- a/include/scopi/vap/vap_fpd.hpp +++ b/include/scopi/vap/vap_fpd.hpp @@ -40,18 +40,7 @@ namespace scopi * @param contacts [in] Array of contacts. */ template - void set_a_priori_velocity_impl(scopi_container& particles, const Contacts& contacts); - - /** - * @brief Constructor. - * - * @param Nactive Number of active particles. - * @param active_ptr Index of the first active particle. - * @param nb_parts Number of objects. - * @param dt Time step. - * @param params Parameters (for compatibility). - */ - vap_fpd(std::size_t Nactive, std::size_t active_ptr, std::size_t nb_parts, double dt); + void set_a_priori_velocity_impl(double dt, scopi_container& particles, const Contacts& contacts); }; /** @@ -78,12 +67,14 @@ namespace scopi type::moment_t<3> cross_product_vap_fpd(const scopi_container<3>& particles, std::size_t i); template - void vap_fpd::set_a_priori_velocity_impl(scopi_container& particles, const Contacts&) + void vap_fpd::set_a_priori_velocity_impl(double dt, scopi_container& particles, const Contacts&) { + auto active_ptr = particles.nb_inactive(); + auto nb_active = particles.nb_active(); #pragma omp parallel for - for (std::size_t i = m_active_ptr; i < m_active_ptr + m_Nactive; ++i) + for (std::size_t i = active_ptr; i < active_ptr + nb_active; ++i) { - particles.v()(i) += m_dt * particles.f()(i) / particles.m()(i); + particles.v()(i) += dt * particles.f()(i) / particles.m()(i); // check cross_product (division by J in the formula missing) and add a torque particles.omega()(i) += cross_product_vap_fpd(particles, i); } diff --git a/include/scopi/vap/vap_projection.hpp b/include/scopi/vap/vap_projection.hpp index 7310bc5d..b3b89f0c 100644 --- a/include/scopi/vap/vap_projection.hpp +++ b/include/scopi/vap/vap_projection.hpp @@ -37,18 +37,7 @@ namespace scopi * @param contacts [in] Array of contacts. */ template - void set_a_priori_velocity_impl(scopi_container& particles, const Contacts& contacts_pos); - - /** - * @brief Constructor. - * - * @param Nactive Number of active particles. - * @param active_ptr Index of the first active particle. - * @param nb_parts Number of objects. - * @param dt Time step. - * @param params Parameters (for compatibility). - */ - vap_projection(std::size_t Nactive, std::size_t active_ptr, std::size_t nb_parts, double dt); + void set_a_priori_velocity_impl(double dt, scopi_container& particles, const Contacts& contacts_pos); /** * @brief Update \c u and \c w. @@ -75,16 +64,19 @@ namespace scopi }; template - void vap_projection::set_a_priori_velocity_impl(scopi_container& particles, const Contacts&) + void vap_projection::set_a_priori_velocity_impl(double, scopi_container& particles, const Contacts&) { + auto active_ptr = particles.nb_inactive(); + auto nb_active = particles.nb_active(); + #pragma omp parallel for - for (std::size_t i = 0; i < this->m_Nactive; ++i) + for (std::size_t i = 0; i < nb_active; ++i) { for (std::size_t d = 0; d < dim; ++d) { - particles.vd()(i + this->m_active_ptr)(d) = m_u(i, d); + particles.vd()(i + active_ptr)(d) = m_u(i, d); } - particles.omega()(i + this->m_active_ptr) = m_w(i, 2); + particles.omega()(i + active_ptr) = m_w(i, 2); } } diff --git a/readme.md b/readme.md index b6320372..71e53185 100644 --- a/readme.md +++ b/readme.md @@ -30,7 +30,7 @@ -The open source software SCoPI written in C++ is designed for the Simulation of Interacting Particle Collections. It enables 2D and 3D simulation of particles (with convex shape), interacting with a surrounding medium (obstacles of various types, moving or not, gravity force, possible coupling with a fluid solver) as well as with other particles (inter-particle forces, contacts). The contact model used is of the Contact Dynamics type [^1] [^2]. It can be dry (inelastic, with or without friction) or viscous (modeling the force of lubrication in a viscous fluid [^3]). The algorithms used are implicit and reduce, at each time step, to the solution of a constrained convex problem. The latter can be solved either using the Mosek library, or via projected gradient or accelerated projected gradient algorithms. In the context of an inelastic contact model, it has been used, for example, to study the local rheology of a granular material in front of a towed sphere [^4]. +The open source software SCoPI written in C++ is designed for the Simulation of Interacting Particle Collections. It enables 2D and 3D simulation of particles (with convex regular shape), interacting with a surrounding medium (obstacles of various types, moving or not, gravity force, possible coupling with a fluid solver) as well as with other particles (inter-particle forces, contacts). The contact model used is of the Contact Dynamics type [^1] [^2]. It can be dry (inelastic, with or without friction) or viscous (modeling the force of lubrication in a viscous fluid [^3]). The algorithms used are implicit and reduce, at each time step, to the solution of a constrained convex problem. The latter can be solved either using the Mosek library, or via projected gradient or accelerated projected gradient algorithms. In the context of an inelastic contact model, it has been used, for example, to study the local rheology of a granular material in front of a towed sphere [^4].
Table of Contents @@ -73,7 +73,7 @@ In this section, we propose an example of the simulation of two spheres with opp We first specify the dimension of our simulation: here is a 2D problem. Then we create two particles represented by spheres. The first argument is the center and the second one is the radius. - SCoPI has other shapes such as superellipsoid, plan, segment and other shapes which are combinations of these. + SCoPI has other shapes such as superellipsoid, plane, segment and other shapes which are combinations of these. - Create the container with all the particles @@ -107,10 +107,10 @@ In this section, we propose an example of the simulation of two spheres with opp ```cpp double dt = 0.005 - scopi::ScopiSolver solver(particles, dt); + scopi::ScopiSolver solver(particles); std::size_t total_it = 100; - solver.run(total_it); + solver.run(dt, total_it); ``` That's it! You have simulated no friction contact of two spheres which is the default behavior. @@ -150,11 +150,11 @@ In this section, we propose an example of the simulation of two spheres with opp double dt = 0.005; - scopi::ScopiSolver solver(particles, dt); + scopi::ScopiSolver solver(particles); SCOPI_PARSE(argc, argv); std::size_t total_it = 100; - solver.run(total_it); + solver.run(dt, total_it); return 0; } @@ -172,6 +172,8 @@ In this section, we propose an example of the simulation of two spheres with opp target_link_librairies(two_spheres scopi) ``` + The files must be in the same directory and outside the core scopi project. + Then, you can compile this C++ script using the bash command lines ```bash @@ -197,7 +199,7 @@ If you want to learn more about scopi skills by looking at examples, we encourag ## Features -- Several objects: sphere, superellipsoid, plan, segment, worms... +- Several objects: sphere, superellipsoid, plane, segment, worms... - Several methods to compute the neighbors of an object: brute force or kd tree - Several kinds of contact: no friction, friction, viscous, friction and viscous - Several methods to solve the problem: projected gradient descent, adaptive projected gradient method diff --git a/src/problems/DryWithFrictionBase.cpp b/src/problems/DryWithFrictionBase.cpp index 81ab3f46..df5db7b7 100644 --- a/src/problems/DryWithFrictionBase.cpp +++ b/src/problems/DryWithFrictionBase.cpp @@ -37,7 +37,7 @@ namespace scopi } std::pair, double> - analytical_solution_sphere_plan_velocity(double alpha, double mu, double t, double r, double g, double y0) + analytical_solution_sphere_plane_velocity(double alpha, double mu, double t, double r, double g, double y0) { double v_normal, omega; double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); diff --git a/src/vap/vap_fixed.cpp b/src/vap/vap_fixed.cpp index 27a7c248..dd51e529 100644 --- a/src/vap/vap_fixed.cpp +++ b/src/vap/vap_fixed.cpp @@ -3,8 +3,5 @@ namespace scopi { - vap_fixed::vap_fixed(std::size_t Nactive, std::size_t active_ptr, std::size_t, double dt) - : base_type(Nactive, active_ptr, dt) - { - } + } diff --git a/src/vap/vap_fpd.cpp b/src/vap/vap_fpd.cpp index 2c6ba04f..0d92e9c3 100644 --- a/src/vap/vap_fpd.cpp +++ b/src/vap/vap_fpd.cpp @@ -3,11 +3,6 @@ namespace scopi { - vap_fpd::vap_fpd(std::size_t Nactive, std::size_t active_ptr, std::size_t, double dt) - : base_type(Nactive, active_ptr, dt) - { - } - type::moment_t<2> cross_product_vap_fpd(const scopi_container<2>&, std::size_t) { return 0.; diff --git a/src/vap/vap_projection.cpp b/src/vap/vap_projection.cpp index 3c31fd42..81f061df 100644 --- a/src/vap/vap_projection.cpp +++ b/src/vap/vap_projection.cpp @@ -3,11 +3,6 @@ namespace scopi { - vap_projection::vap_projection(std::size_t Nactive, std::size_t active_ptr, std::size_t, double dt) - : base_type(Nactive, active_ptr, dt) - { - } - void vap_projection::set_u_w(const xt::xtensor& u, const xt::xtensor& w) { m_u = u; diff --git a/test/analytical_solution.cpp b/test/analytical_solution.cpp index 1d08a511..59e8da25 100644 --- a/test/analytical_solution.cpp +++ b/test/analytical_solution.cpp @@ -4,7 +4,7 @@ namespace scopi { /// Dry without friction std::pair, double> - analytical_solution_sphere_plan_velocity_no_friction(double alpha, double t, double r, double g, double y0) + analytical_solution_sphere_plane_velocity_no_friction(double alpha, double t, double r, double g, double y0) { const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); type::position_t<2> v; @@ -20,7 +20,7 @@ namespace scopi return std::make_pair(v, 0.); } - std::pair, double> analytical_solution_sphere_plan_no_friction(double alpha, double t, double r, double g, double y0) + std::pair, double> analytical_solution_sphere_plane_no_friction(double alpha, double t, double r, double g, double y0) { const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); type::position_t<2> x; @@ -38,7 +38,7 @@ namespace scopi /// Visous without friction std::pair, double> - analytical_solution_sphere_plan_velocity_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv) + analytical_solution_sphere_plane_velocity_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv) { type::position_t<2> v; const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); @@ -80,7 +80,7 @@ namespace scopi } std::pair, double> - analytical_solution_sphere_plan_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv) + analytical_solution_sphere_plane_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv) { type::position_t<2> x; const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); @@ -128,7 +128,7 @@ namespace scopi /// Friction std::pair, double> - analytical_solution_sphere_plan_friction(double alpha, double mu, double t, double r, double g, double y0) + analytical_solution_sphere_plane_friction(double alpha, double mu, double t, double r, double g, double y0) { const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); type::position_t<2> x; @@ -160,7 +160,7 @@ namespace scopi } std::pair, double> - analytical_solution_sphere_plan_velocity_friction(double alpha, double mu, double t, double r, double g, double y0) + analytical_solution_sphere_plane_velocity_friction(double alpha, double mu, double t, double r, double g, double y0) { const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); type::position_t<2> v; @@ -191,14 +191,14 @@ namespace scopi } /// Viscous with friction - std::pair, double> analytical_solution_sphere_plan_velocity_viscous_friction(double alpha, - double mu, - double t, - double r, - double g, - double y0, - double gamma_min, - double t_inv) + std::pair, double> analytical_solution_sphere_plane_velocity_viscous_friction(double alpha, + double mu, + double t, + double r, + double g, + double y0, + double gamma_min, + double t_inv) { double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); double t_gamma_min = -gamma_min / (g * std::cos(alpha)); @@ -295,7 +295,7 @@ namespace scopi } std::pair, double> - analytical_solution_sphere_plan_viscous_friction(double alpha, double mu, double t, double r, double g, double y0, double gamma_min, double t_inv) + analytical_solution_sphere_plane_viscous_friction(double alpha, double mu, double t, double r, double g, double y0, double gamma_min, double t_inv) { double angle; double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); diff --git a/test/analytical_solution.hpp b/test/analytical_solution.hpp index 0306269b..bdf1efdb 100644 --- a/test/analytical_solution.hpp +++ b/test/analytical_solution.hpp @@ -9,33 +9,34 @@ namespace scopi { std::pair, double> - analytical_solution_sphere_plan_velocity_no_friction(double alpha, double t, double r, double g, double y0); - std::pair, double> analytical_solution_sphere_plan_no_friction(double alpha, double t, double r, double g, double y0); + analytical_solution_sphere_plane_velocity_no_friction(double alpha, double t, double r, double g, double y0); + std::pair, double> + analytical_solution_sphere_plane_no_friction(double alpha, double t, double r, double g, double y0); std::pair, double> - analytical_solution_sphere_plan_velocity_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv); + analytical_solution_sphere_plane_velocity_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv); std::pair, double> - analytical_solution_sphere_plan_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv); + analytical_solution_sphere_plane_viscous(double alpha, double t, double r, double g, double y0, double gamma_min, double t_inv); std::pair, double> - analytical_solution_sphere_plan_friction(double alpha, double mu, double t, double r, double g, double y0); + analytical_solution_sphere_plane_friction(double alpha, double mu, double t, double r, double g, double y0); std::pair, double> - analytical_solution_sphere_plan_velocity_friction(double alpha, double mu, double t, double r, double g, double y0); + analytical_solution_sphere_plane_velocity_friction(double alpha, double mu, double t, double r, double g, double y0); - std::pair, double> analytical_solution_sphere_plan_velocity_viscous_friction(double alpha, - double mu, - double t, - double r, - double g, - double y0, - double gamma_min, - double t_inv); - std::pair, double> analytical_solution_sphere_plan_viscous_friction(double alpha, - double mu, - double t, - double r, - double g, - double y0, - double gamma_min, - double t_inv); + std::pair, double> analytical_solution_sphere_plane_velocity_viscous_friction(double alpha, + double mu, + double t, + double r, + double g, + double y0, + double gamma_min, + double t_inv); + std::pair, double> analytical_solution_sphere_plane_viscous_friction(double alpha, + double mu, + double t, + double r, + double g, + double y0, + double gamma_min, + double t_inv); } diff --git a/test/references/sphere_plan_friction_mu01.json b/test/references/sphere_plane_friction_mu01.json similarity index 98% rename from test/references/sphere_plan_friction_mu01.json rename to test/references/sphere_plane_friction_mu01.json index 4d79ee21..8c5b0c8b 100644 --- a/test/references/sphere_plan_friction_mu01.json +++ b/test/references/sphere_plane_friction_mu01.json @@ -43,7 +43,7 @@ 0.8660254037844386, 0.5 ], - "type": "plan" + "type": "plane" }, { "id": 1, diff --git a/test/references/sphere_plan_viscosity.json b/test/references/sphere_plane_viscosity.json similarity index 98% rename from test/references/sphere_plan_viscosity.json rename to test/references/sphere_plane_viscosity.json index 0c601984..621756c2 100644 --- a/test/references/sphere_plan_viscosity.json +++ b/test/references/sphere_plane_viscosity.json @@ -38,7 +38,7 @@ 0.7071067811865476, 0.7071067811865475 ], - "type": "plan" + "type": "plane" }, { "id": 1, diff --git a/test/references/sphere_plan_viscosity_friction.json b/test/references/sphere_plane_viscosity_friction.json similarity index 98% rename from test/references/sphere_plan_viscosity_friction.json rename to test/references/sphere_plane_viscosity_friction.json index 10db39f2..c2164b35 100644 --- a/test/references/sphere_plan_viscosity_friction.json +++ b/test/references/sphere_plane_viscosity_friction.json @@ -38,7 +38,7 @@ 1.0, 2.220446049250313e-16 ], - "type": "plan" + "type": "plane" }, { "id": 1, diff --git a/test/references/sphere_plan_viscosity_friction_vertical.json b/test/references/sphere_plane_viscosity_friction_vertical.json similarity index 97% rename from test/references/sphere_plan_viscosity_friction_vertical.json rename to test/references/sphere_plane_viscosity_friction_vertical.json index 06422231..5557f15c 100644 --- a/test/references/sphere_plan_viscosity_friction_vertical.json +++ b/test/references/sphere_plane_viscosity_friction_vertical.json @@ -38,7 +38,7 @@ 1.0, 2.220446049250313e-16 ], - "type": "plan" + "type": "plane" }, { "id": 1, diff --git a/test/test_closest_points.cpp b/test/test_closest_points.cpp index 7cf809da..91a2ba83 100644 --- a/test/test_closest_points.cpp +++ b/test/test_closest_points.cpp @@ -250,8 +250,8 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.4)); } - // distance sphere - plan - TEST_CASE("sphere_plan_2d") + // distance sphere - plane + TEST_CASE("sphere_plane_2d") { constexpr std::size_t dim = 2; sphere s( @@ -259,7 +259,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0} }, @@ -276,7 +276,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_2d_rotation_30_deg") + TEST_CASE("sphere_plane_2d_rotation_30_deg") { constexpr std::size_t dim = 2; sphere s( @@ -287,7 +287,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot} }, @@ -304,7 +304,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_2d_rotation_90_deg") + TEST_CASE("sphere_plane_2d_rotation_90_deg") { constexpr std::size_t dim = 2; sphere s( @@ -312,7 +312,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0., -0.2} }, @@ -329,7 +329,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - TEST_CASE("sphere_plan_2d_dispatch") + TEST_CASE("sphere_plane_2d_dispatch") { constexpr std::size_t dim = 2; sphere s( @@ -337,7 +337,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0} }, @@ -358,7 +358,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_2d_dispatch_rotation_30_deg") + TEST_CASE("sphere_plane_2d_dispatch_rotation_30_deg") { constexpr std::size_t dim = 2; sphere s( @@ -369,7 +369,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot} }, @@ -390,7 +390,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_2d_dispatch_rotation_90_deg") + TEST_CASE("sphere_plane_2d_dispatch_rotation_90_deg") { constexpr std::size_t dim = 2; sphere s( @@ -398,7 +398,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0., -0.2} }, @@ -419,7 +419,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - TEST_CASE("sphere_plan_3d") + TEST_CASE("sphere_plane_3d") { constexpr std::size_t dim = 3; sphere s( @@ -427,7 +427,7 @@ namespace scopi {0.0, 0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0, 0.0} }, @@ -447,7 +447,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_3d_rotation_30_deg") + TEST_CASE("sphere_plane_3d_rotation_30_deg") { constexpr std::size_t dim = 3; sphere s( @@ -458,7 +458,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot, 0.} }, @@ -478,7 +478,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_3d_rotation_90_deg") + TEST_CASE("sphere_plane_3d_rotation_90_deg") { constexpr std::size_t dim = 3; sphere s( @@ -488,7 +488,7 @@ namespace scopi 0, } }, 0.1); - plan p( + plane p( { {0., -0.2, 0.} }, @@ -508,7 +508,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - TEST_CASE("sphere_plan_3d_dispatch") + TEST_CASE("sphere_plane_3d_dispatch") { constexpr std::size_t dim = 3; sphere s( @@ -516,7 +516,7 @@ namespace scopi {0.0, 0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0, 0.0} }, @@ -540,7 +540,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_3d_dispatch_rotation_30_deg") + TEST_CASE("sphere_plane_3d_dispatch_rotation_30_deg") { constexpr std::size_t dim = 3; sphere s( @@ -551,7 +551,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot, 0.} }, @@ -575,7 +575,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("sphere_plan_3d_dispatch_rotation_90_deg") + TEST_CASE("sphere_plane_3d_dispatch_rotation_90_deg") { constexpr std::size_t dim = 3; sphere s( @@ -585,7 +585,7 @@ namespace scopi 0, } }, 0.1); - plan p( + plane p( { {0., -0.2, 0.} }, @@ -609,8 +609,8 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - // distance plan - sphere - TEST_CASE("plan_sphere_2d") + // distance plane - sphere + TEST_CASE("plane_sphere_2d") { constexpr std::size_t dim = 2; sphere s( @@ -618,7 +618,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0} }, @@ -635,7 +635,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_2d_rotation_30_deg") + TEST_CASE("plane_sphere_2d_rotation_30_deg") { constexpr std::size_t dim = 2; sphere s( @@ -646,7 +646,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot} }, @@ -663,7 +663,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_2d_rotation_90_deg") + TEST_CASE("plane_sphere_2d_rotation_90_deg") { constexpr std::size_t dim = 2; sphere s( @@ -671,7 +671,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0., -0.2} }, @@ -688,7 +688,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - TEST_CASE("plan_sphere_2d_dispatch") + TEST_CASE("plane_sphere_2d_dispatch") { constexpr std::size_t dim = 2; sphere s( @@ -696,7 +696,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0} }, @@ -717,7 +717,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_2d_dispatch_rotation_30_deg") + TEST_CASE("plane_sphere_2d_dispatch_rotation_30_deg") { constexpr std::size_t dim = 2; sphere s( @@ -728,7 +728,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot} }, @@ -749,7 +749,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_2d_dispatch_rotation_90_deg") + TEST_CASE("plane_sphere_2d_dispatch_rotation_90_deg") { constexpr std::size_t dim = 2; sphere s( @@ -757,7 +757,7 @@ namespace scopi {0.0, 0.0} }, 0.1); - plan p( + plane p( { {0., -0.2} }, @@ -778,7 +778,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - TEST_CASE("plan_sphere_3d") + TEST_CASE("plane_sphere_3d") { constexpr std::size_t dim = 3; sphere s( @@ -786,7 +786,7 @@ namespace scopi {0.0, 0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0, 0.0} }, @@ -806,7 +806,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_3d_rotation_30_deg") + TEST_CASE("plane_sphere_3d_rotation_30_deg") { constexpr std::size_t dim = 3; sphere s( @@ -817,7 +817,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot, 0.} }, @@ -837,7 +837,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_3d_rotation_90_deg") + TEST_CASE("plane_sphere_3d_rotation_90_deg") { constexpr std::size_t dim = 3; sphere s( @@ -847,7 +847,7 @@ namespace scopi 0, } }, 0.1); - plan p( + plane p( { {0., -0.2, 0.} }, @@ -867,7 +867,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.1)); } - TEST_CASE("plan_sphere_3d_dispatch") + TEST_CASE("plane_sphere_3d_dispatch") { constexpr std::size_t dim = 3; sphere s( @@ -875,7 +875,7 @@ namespace scopi {0.0, 0.0, 0.0} }, 0.1); - plan p( + plane p( { {0.3, 0.0, 0.0} }, @@ -899,7 +899,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_3d_dispatch_rotation_30_deg") + TEST_CASE("plane_sphere_3d_dispatch_rotation_30_deg") { constexpr std::size_t dim = 3; sphere s( @@ -910,7 +910,7 @@ namespace scopi double dist = 0.3; double cosRot = std::sqrt(3.) / 2.; double sinRot = 1. / 2.; - plan p( + plane p( { {dist * cosRot, dist * sinRot, 0.} }, @@ -934,7 +934,7 @@ namespace scopi REQUIRE(out.dij == doctest::Approx(0.2)); } - TEST_CASE("plan_sphere_3d_dispatch_rotation_90_deg") + TEST_CASE("plane_sphere_3d_dispatch_rotation_90_deg") { constexpr std::size_t dim = 3; sphere s( @@ -944,7 +944,7 @@ namespace scopi 0, } }, 0.1); - plan p( + plane p( { {0., -0.2, 0.} }, @@ -1771,8 +1771,8 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.5)); // } - // // distance superellipsoid - plan - // TEST_CASE("superellipsoid_plan_2d" * doctest::may_fail(true)) + // // distance superellipsoid - plane + // TEST_CASE("superellipsoid_plane_2d" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -1783,7 +1783,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2}}, // 1); - // plan p( + // plane p( // { // {0.3, 0.0} // }, @@ -1800,7 +1800,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_2d_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_2d_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -1814,7 +1814,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot} // }, @@ -1831,7 +1831,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_2d_dispatch" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_2d_dispatch" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -1842,7 +1842,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2}}, // 1); - // plan p( + // plane p( // { // {0.3, 0.0} // }, @@ -1863,7 +1863,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_2d_dispatch_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_2d_dispatch_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -1877,7 +1877,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot} // }, @@ -1898,7 +1898,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_3d" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_3d" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -1909,7 +1909,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2, 0.3}}, // {1, 1}); - // plan p( + // plane p( // { // {0.3, 0.0, 0.0} // }, @@ -1929,7 +1929,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_3d_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_3d_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -1943,7 +1943,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot, 0.} // }, @@ -1963,7 +1963,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_3d_dispatch" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_3d_dispatch" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -1974,7 +1974,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2, 0.3}}, // {1, 1}); - // plan p( + // plane p( // { // {0.3, 0.0, 0.0} // }, @@ -1998,7 +1998,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("superellipsoid_plan_3d_dispatch_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("superellipsoid_plane_3d_dispatch_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -2012,7 +2012,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot, 0.} // }, @@ -2036,8 +2036,8 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // // distance plan - superellipsoid - // TEST_CASE("plan_superellipsoid_2d" * doctest::may_fail(true)) + // // distance plane - superellipsoid + // TEST_CASE("plane_superellipsoid_2d" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -2048,7 +2048,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2}}, // 1); - // plan p( + // plane p( // { // {0.3, 0.0} // }, @@ -2065,7 +2065,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_2d_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_2d_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -2079,7 +2079,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot} // }, @@ -2096,7 +2096,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_2d_dispatch" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_2d_dispatch" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -2107,7 +2107,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2}}, // 1); - // plan p( + // plane p( // { // {0.3, 0.0} // }, @@ -2128,7 +2128,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_2d_dispatch_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_2d_dispatch_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 2; @@ -2142,7 +2142,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot} // }, @@ -2163,7 +2163,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_3d" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_3d" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -2174,7 +2174,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2, 0.3}}, // {1, 1}); - // plan p( + // plane p( // { // {0.3, 0.0, 0.0} // }, @@ -2194,7 +2194,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_3d_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_3d_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -2208,7 +2208,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot, 0.} // }, @@ -2228,7 +2228,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_3d_dispatch" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_3d_dispatch" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -2239,7 +2239,7 @@ namespace scopi // {quaternion(0.)}, // {{0.1, 0.2, 0.3}}, // {1, 1}); - // plan p( + // plane p( // { // {0.3, 0.0, 0.0} // }, @@ -2263,7 +2263,7 @@ namespace scopi // REQUIRE(out.dij == doctest::Approx(0.2)); // } - // TEST_CASE("plan_superellipsoid_3d_dispatch_rotation_30_deg" * doctest::may_fail(true)) + // TEST_CASE("plane_superellipsoid_3d_dispatch_rotation_30_deg" * doctest::may_fail(true)) // { // // FIXME Newton does not converge // constexpr std::size_t dim = 3; @@ -2277,7 +2277,7 @@ namespace scopi // double dist = 0.3; // double cosRot = std::sqrt(3.) / 2.; // double sinRot = 1. / 2.; - // plan p( + // plane p( // { // {dist * cosRot, dist * sinRot, 0.} // }, diff --git a/test/test_friction.cpp b/test/test_friction.cpp index 7f193c48..936e7b57 100644 --- a/test/test_friction.cpp +++ b/test/test_friction.cpp @@ -43,11 +43,11 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.default_contact_property.mu = mu; // params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_spheres_asymmetrical_friction.json", tolerance)); } @@ -102,16 +102,16 @@ namespace scopi } } - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.default_contact_property.mu = mu; params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0099.json", "../test/references/2d_case_spheres_friction.json", tolerance)); } - TEST_CASE_TEMPLATE_DEFINE("sphere inclined plan friction", SolverType, sphere_inclined_plan_friction) + TEST_CASE_TEMPLATE_DEFINE("sphere inclined plane friction", SolverType, sphere_inclined_plane_friction) { std::tuple data; std::vector> data_container( @@ -136,7 +136,7 @@ namespace scopi double alpha = std::get<1>(data); auto prop = property().mass(1.).moment_inertia(1. * radius * radius / 2.); - plan p( + plane p( { {0., 0.} }, @@ -154,11 +154,11 @@ namespace scopi {0., -g} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.default_contact_property.mu = mu; params.solver_params.output_frequency = std::size_t(-1); - solver.run(total_it); + solver.run(dt, total_it); auto pos = particles.pos(); auto q = particles.q(); @@ -169,7 +169,7 @@ namespace scopi double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); auto v = particles.v(); auto omega = particles.omega(); - tmp = analytical_solution_sphere_plan_velocity(alpha, mu, dt * (total_it + 1), radius, g, h); + tmp = analytical_solution_sphere_plane_velocity(alpha, mu, dt * (total_it + 1), radius, g, h); auto v_analytical = tmp.first; double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); auto omega_analytical = tmp.second; @@ -183,6 +183,6 @@ namespace scopi TEST_CASE_TEMPLATE_APPLY(two_spheres_asymetrical_friction, solver_dry_with_friction_t<2, vap_fixed>); TEST_CASE_TEMPLATE_APPLY(critical_2d_spheres_friction, solver_dry_with_friction_t<2, vap_fixed>); - TEST_CASE_TEMPLATE_APPLY(sphere_inclined_plan_friction, solver_dry_with_friction_t<2, vap_fpd>); + TEST_CASE_TEMPLATE_APPLY(sphere_inclined_plane_friction, solver_dry_with_friction_t<2, vap_fpd>); } diff --git a/test/test_gradient.cpp b/test/test_gradient.cpp index f9458ff6..b1e8b12f 100644 --- a/test/test_gradient.cpp +++ b/test/test_gradient.cpp @@ -1,7 +1,7 @@ #include #include -#include +#include #include #include #include @@ -19,8 +19,8 @@ namespace scopi double tol = 1e-5; double tol_analytical = 1e-12; - /// TESTS SPHERES - PLAN - TEST_CASE("sphere - plan without friction") + /// TESTS SPHERES - PLANE + TEST_CASE("sphere - plane without friction") { constexpr std::size_t dim = 2; @@ -32,7 +32,7 @@ namespace scopi std::size_t total_it = 1000; double h = 2. * radius; double alpha = PI / 4; - plan plan( + plane plane( { {0., 0.} }, @@ -42,7 +42,7 @@ namespace scopi {0, h} }, radius); - particles.push_back(plan, property().deactivate()); + particles.push_back(plane, property().deactivate()); particles.push_back(sphere, property() .mass(1) @@ -54,7 +54,7 @@ namespace scopi using problem_t = NoFriction; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); params.optim_params.tolerance = 1e-7; @@ -63,19 +63,19 @@ namespace scopi params.optim_params.dynamic_descent = true; params.solver_params.output_frequency = total_it - 1; - params.solver_params.filename = "sphere_plan_nofriction"; - solver.run(total_it); + params.solver_params.filename = "sphere_plane_nofriction"; + solver.run(dt, total_it); auto pos = particles.pos(); auto q = particles.q(); - auto tmp = analytical_solution_sphere_plan_no_friction(alpha, dt * (total_it), radius, g, h); + auto tmp = analytical_solution_sphere_plane_no_friction(alpha, dt * (total_it), radius, g, h); auto pos_analytical = tmp.first; auto q_analytical = quaternion(tmp.second); double error_pos = xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical); double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); auto v = particles.v(); auto omega = particles.omega(); - tmp = analytical_solution_sphere_plan_velocity_no_friction(alpha, dt * (total_it), radius, g, h); + tmp = analytical_solution_sphere_plane_velocity_no_friction(alpha, dt * (total_it), radius, g, h); auto v_analytical = tmp.first; double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); @@ -83,10 +83,10 @@ namespace scopi REQUIRE(error_q <= doctest::Approx(1e-3)); REQUIRE(error_v <= doctest::Approx(tol_analytical)); REQUIRE(omega(1) == doctest::Approx(0.)); - // CHECK(diffFile("./Results/sphere_plan_nofriction_0002.json", "../test/references/sphere_plan_nofriction.json", tol)); + // CHECK(diffFile("./Results/sphere_plane_nofriction_0002.json", "../test/references/sphere_plane_nofriction.json", tol)); } - TEST_CASE("sphere - plan viscous without friction") + TEST_CASE("sphere - plane viscous without friction") { constexpr std::size_t dim = 2; @@ -98,7 +98,7 @@ namespace scopi std::size_t total_it = 500; double h = 2. * radius; double alpha = PI / 4; - plan plan( + plane plane( { {0., 0.} }, @@ -108,7 +108,7 @@ namespace scopi {0, h} }, radius); - particles.push_back(plan, property().deactivate()); + particles.push_back(plane, property().deactivate()); particles.push_back(sphere, property() .mass(1) @@ -120,7 +120,7 @@ namespace scopi using problem_t = Viscous; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); params.optim_params.tolerance = 1e-7; @@ -129,27 +129,27 @@ namespace scopi params.optim_params.dynamic_descent = true; params.solver_params.output_frequency = 1; - params.solver_params.filename = "sphere_plan_viscous"; + params.solver_params.filename = "sphere_plane_viscous"; // double gamma = 0; double gamma_min = -2.; // double gamma_tol = 1e-6; - solver.run(total_it); + solver.run(dt, total_it); particles.f()(1) = {-g * sin(alpha), g * cos(alpha)}; double total_it_2 = 500; - solver.run(total_it_2); + solver.run(dt, total_it_2); auto pos = particles.pos(); auto q = particles.q(); - auto tmp = analytical_solution_sphere_plan_viscous(alpha, dt * (total_it + total_it_2), radius, g, h, gamma_min, dt * (total_it)); + auto tmp = analytical_solution_sphere_plane_viscous(alpha, dt * (total_it + total_it_2), radius, g, h, gamma_min, dt * (total_it)); auto pos_analytical = tmp.first; auto q_analytical = quaternion(tmp.second); double error_pos = xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical); double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); auto v = particles.v(); auto omega = particles.omega(); - tmp = analytical_solution_sphere_plan_velocity_viscous(alpha, dt * (total_it + total_it_2), radius, g, h, gamma_min, dt * (total_it)); + tmp = analytical_solution_sphere_plane_velocity_viscous(alpha, dt * (total_it + total_it_2), radius, g, h, gamma_min, dt * (total_it)); auto v_analytical = tmp.first; double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); @@ -157,16 +157,16 @@ namespace scopi REQUIRE(error_q <= doctest::Approx(1e-3)); REQUIRE(error_v <= doctest::Approx(1e-3)); REQUIRE(omega(1) == doctest::Approx(0.)); - // CHECK(diffFile("./Results/sphere_plan_viscous_0005.json", "../test/references/sphere_plan_viscous.json", tol)); + // CHECK(diffFile("./Results/sphere_plane_viscous_0005.json", "../test/references/sphere_plane_viscous.json", tol)); } - TEST_CASE("sphere - plan friction no fixed point") + TEST_CASE("sphere - plane friction no fixed point") { constexpr std::size_t dim = 2; scopi_container particles; - plan plan( + plane plane( { {0., 0.} }, @@ -180,7 +180,7 @@ namespace scopi }, rr); - particles.push_back(plan, property().deactivate()); + particles.push_back(plane, property().deactivate()); particles.push_back(sphere, property() .mass(1) @@ -195,7 +195,7 @@ namespace scopi using problem_t = Friction; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); @@ -205,16 +205,16 @@ namespace scopi params.optim_params.dynamic_descent = true; params.solver_params.output_frequency = Tf / dt - 1; - params.solver_params.filename = "sphere_plan_friction_mu01"; + params.solver_params.filename = "sphere_plane_friction_mu01"; // mu01 params.default_contact_property.mu = 0.1; - solver.run(Tf / dt); + solver.run(dt, Tf / dt); - CHECK(diffFile("./Results/sphere_plan_friction_mu01_0200.json", "../test/references/sphere_plan_friction_mu01.json", tol)); + CHECK(diffFile("./Results/sphere_plane_friction_mu01_0200.json", "../test/references/sphere_plane_friction_mu01.json", tol)); } - TEST_CASE("sphere - plan friction fixed point") + TEST_CASE("sphere - plane friction fixed point") { constexpr std::size_t dim = 2; @@ -226,7 +226,7 @@ namespace scopi std::size_t total_it = 500; double h = 2. * radius; double alpha = PI / 4; - plan plan( + plane plane( { {0., 0.} }, @@ -236,7 +236,7 @@ namespace scopi {0, h} }, radius); - particles.push_back(plan, property().deactivate()); + particles.push_back(plane, property().deactivate()); particles.push_back(sphere, property() .mass(1) @@ -248,7 +248,7 @@ namespace scopi using problem_t = FrictionFixedPoint; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); @@ -259,22 +259,22 @@ namespace scopi double mu = 0.5; params.solver_params.output_frequency = total_it - 1; - params.solver_params.filename = "sphere_plan_friction_fixed_point"; + params.solver_params.filename = "sphere_plane_friction_fixed_point"; params.default_contact_property.mu = mu; params.default_contact_property.fixed_point_tol = 1e-6; params.default_contact_property.fixed_point_max_iter = 1000; - solver.run(total_it); + solver.run(dt, total_it); auto pos = particles.pos(); auto q = particles.q(); - auto tmp = analytical_solution_sphere_plan_friction(alpha, mu, dt * (total_it), radius, g, h); + auto tmp = analytical_solution_sphere_plane_friction(alpha, mu, dt * (total_it), radius, g, h); auto pos_analytical = tmp.first; auto q_analytical = quaternion(tmp.second); double error_pos = xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical); double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); auto v = particles.v(); auto omega = particles.omega(); - tmp = analytical_solution_sphere_plan_velocity_friction(alpha, mu, dt * (total_it), radius, g, h); + tmp = analytical_solution_sphere_plane_velocity_friction(alpha, mu, dt * (total_it), radius, g, h); auto v_analytical = tmp.first; auto omega_analytical = tmp.second; double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); @@ -283,10 +283,11 @@ namespace scopi REQUIRE(error_q <= doctest::Approx(1e-2)); REQUIRE(error_v <= doctest::Approx(tol_analytical)); REQUIRE(error_omega <= doctest::Approx(tol_analytical)); - // CHECK(diffFile("./Results/sphere_plan_friction_fp_mu05_0002.json", "../test/references/sphere_plan_friction_fp_mu05.json", tol)); + // CHECK(diffFile("./Results/sphere_plane_friction_fp_mu05_0002.json", "../test/references/sphere_plane_friction_fp_mu05.json", + // tol)); } - TEST_CASE("sphere - plan viscous friction") + TEST_CASE("sphere - plane viscous friction") { constexpr std::size_t dim = 2; @@ -298,7 +299,7 @@ namespace scopi std::size_t total_it = 500; double h = 2. * radius; double alpha = PI / 4; - plan plan( + plane plane( { {0., 0.} }, @@ -308,7 +309,7 @@ namespace scopi {0, h} }, radius); - particles.push_back(plan, property().deactivate()); + particles.push_back(plane, property().deactivate()); particles.push_back(sphere, property() .mass(1) @@ -320,7 +321,7 @@ namespace scopi using problem_t = ViscousFriction; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); params.optim_params.tolerance = 1e-7; @@ -329,7 +330,7 @@ namespace scopi params.optim_params.dynamic_descent = true; params.solver_params.output_frequency = total_it - 1; - params.solver_params.filename = "sphere_plan_viscous_friction_mu05"; + params.solver_params.filename = "sphere_plane_viscous_friction_mu05"; double mu = 0.5; params.default_contact_property.mu = mu; @@ -340,35 +341,35 @@ namespace scopi params.default_contact_property.gamma_min = gamma_min; params.default_contact_property.gamma_tol = 1e-6; - solver.run(total_it); + solver.run(dt, total_it); particles.f()(1) = {-g * sin(alpha), g * cos(alpha)}; double total_it_2 = 100; - solver.run(total_it_2); + solver.run(dt, total_it_2); auto pos = particles.pos(); auto q = particles.q(); - auto tmp = analytical_solution_sphere_plan_viscous_friction(alpha, - mu, - dt * (total_it + total_it_2), - radius, - g, - h, - gamma_min, - dt * (total_it)); + auto tmp = analytical_solution_sphere_plane_viscous_friction(alpha, + mu, + dt * (total_it + total_it_2), + radius, + g, + h, + gamma_min, + dt * (total_it)); auto pos_analytical = tmp.first; auto q_analytical = quaternion(tmp.second); double error_pos = xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical); double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); auto v = particles.v(); auto omega = particles.omega(); - tmp = analytical_solution_sphere_plan_velocity_viscous_friction(alpha, - mu, - dt * (total_it + total_it_2), - radius, - g, - h, - gamma_min, - dt * (total_it)); + tmp = analytical_solution_sphere_plane_velocity_viscous_friction(alpha, + mu, + dt * (total_it + total_it_2), + radius, + g, + h, + gamma_min, + dt * (total_it)); auto v_analytical = tmp.first; auto omega_analytical = tmp.second; double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); @@ -379,8 +380,8 @@ namespace scopi REQUIRE(error_v <= doctest::Approx(1e-3)); REQUIRE(error_omega <= doctest::Approx(1e-3)); - // CHECK(diffFile("./Results/sphere_plan_viscous_friction_mu05_0005.json", - // "../test/references/sphere_plan_viscous_friction_mu05.json", tol)); + // CHECK(diffFile("./Results/sphere_plane_viscous_friction_mu05_0005.json", + // "../test/references/sphere_plane_viscous_friction_mu05.json", tol)); } /// TESTS 3 SPHERES @@ -428,7 +429,7 @@ namespace scopi using problem_t = NoFriction; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); @@ -439,7 +440,7 @@ namespace scopi params.solver_params.output_frequency = Tf / dt - 1; params.solver_params.filename = "3spheres_nofriction"; - solver.run(Tf / dt); + solver.run(dt, Tf / dt); CHECK(diffFile("./Results/3spheres_nofriction_0060.json", "../test/references/3spheres_nofriction.json", tol)); } @@ -486,7 +487,7 @@ namespace scopi using problem_t = Viscous; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); @@ -501,7 +502,7 @@ namespace scopi // double gamma = 0; // double gamma_min = -2.; // double gamma_tol = 1e-6; - solver.run(Tf / dt); + solver.run(dt, Tf / dt); CHECK(diffFile("./Results/3spheres_viscous_0070.json", "../test/references/3spheres_viscous.json", tol)); } @@ -548,7 +549,7 @@ namespace scopi using problem_t = FrictionFixedPoint; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); @@ -564,7 +565,7 @@ namespace scopi params.default_contact_property.mu = 0.5; params.default_contact_property.fixed_point_tol = 1e-6; params.default_contact_property.fixed_point_max_iter = 1000; - solver.run(Tf / dt); + solver.run(dt, Tf / dt); CHECK(diffFile("./Results/3spheres_friction_fp_0070.json", "../test/references/3spheres_friction_fp.json", tol)); } @@ -611,7 +612,7 @@ namespace scopi using problem_t = ViscousFriction; using optim_solver = OptimGradient; using vap_t = vap_fpd; - ScopiSolver solver(particles, dt); + ScopiSolver solver(particles); auto params = solver.get_params(); @@ -629,7 +630,7 @@ namespace scopi params.default_contact_property.gamma = 0; params.default_contact_property.gamma_min = -1.4; params.default_contact_property.gamma_tol = 1e-6; - solver.run(Tf / dt); + solver.run(dt, Tf / dt); CHECK(diffFile("./Results/3spheres_viscous_friction_0070.json", "../test/references/3spheres_viscous_friction.json", tol)); } } diff --git a/test/test_matrices.cpp b/test/test_matrices.cpp index f151ae26..217cfe59 100644 --- a/test/test_matrices.cpp +++ b/test/test_matrices.cpp @@ -56,13 +56,13 @@ namespace scopi {2., 1.} }, 0.5); - scopi::plan plan( + scopi::plane plane( { {0., 0.} }, PI / 2 - PI / 4); - particles.push_back(plan, scopi::property().deactivate()); + particles.push_back(plane, scopi::property().deactivate()); particles.push_back(sphere, scopi::property().mass(1).moment_inertia(0.1)); ContactsParams> params; diff --git a/test/test_obstacles.cpp b/test/test_obstacles.cpp index 668a65a1..9ce57f9a 100644 --- a/test/test_obstacles.cpp +++ b/test/test_obstacles.cpp @@ -7,7 +7,7 @@ #include "utils.hpp" #include -#include +#include #include #include #include @@ -40,7 +40,7 @@ namespace scopi params.solver_params.output_frequency = total_it - 1; } - TEST_CASE_TEMPLATE_DEFINE("sphere plan", SolverType, sphere_plan) + TEST_CASE_TEMPLATE_DEFINE("sphere plane", SolverType, sphere_plane) { static constexpr std::size_t dim = 2; @@ -53,7 +53,7 @@ namespace scopi {0., radius} }, radius); - plan p( + plane p( { {0., 0.} }, @@ -66,8 +66,8 @@ namespace scopi SUBCASE("fixed") { particles.push_back(s, prop); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); check_result_sphere_plan(particles); } @@ -78,14 +78,14 @@ namespace scopi prop.desired_velocity({ {0., -1.} })); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); check_result_sphere_plan(particles); } } - TEST_CASE_TEMPLATE_DEFINE("sphere plan force", SolverType, sphere_plan_force) + TEST_CASE_TEMPLATE_DEFINE("sphere plane force", SolverType, sphere_plane_force) { static constexpr std::size_t dim = 2; @@ -98,7 +98,7 @@ namespace scopi {0., radius} }, radius); - plan p( + plane p( { {0., 0.} }, @@ -112,8 +112,8 @@ namespace scopi prop.force({ {0., -1.} })); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); check_result_sphere_plan(particles); } @@ -162,8 +162,8 @@ namespace scopi SUBCASE("fixed") { particles.push_back(sphere, prop); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); check_result_sphere_sphere(particles); } @@ -174,8 +174,8 @@ namespace scopi prop.desired_velocity({ {0., -1.} })); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); check_result_sphere_sphere(particles); } @@ -207,8 +207,8 @@ namespace scopi prop.force({ {0., -1.} })); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); check_result_sphere_sphere(particles); } @@ -239,16 +239,16 @@ namespace scopi {0., -10.} })); - SolverType solver(particles, dt); + SolverType solver(particles); set_params(solver, total_it); auto params = solver.get_params(); params.solver_params.output_frequency = 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0099.json", "../test/references/obstacles_sphere_sphere_moving.json", tolerance)); } - TEST_CASE_TEMPLATE_DEFINE("sphere inclined plan", SolverType, sphere_inclined_plan) + TEST_CASE_TEMPLATE_DEFINE("sphere inclined plane", SolverType, sphere_inclined_plane) { std::tuple data; std::vector> data_container({std::make_tuple(PI / 6., 0.000998789, 0., 0.0010002), @@ -266,7 +266,7 @@ namespace scopi double alpha = std::get<0>(data); auto prop = property().mass(1.).moment_inertia(1. * radius * radius / 2.); - plan p( + plane p( { {0., 0.} }, @@ -284,20 +284,20 @@ namespace scopi {0., -g} })); - SolverType solver(particles, dt); - solver.run(total_it); + SolverType solver(particles); + solver.run(dt, total_it); set_params(solver, total_it); auto pos = particles.pos(); auto q = particles.q(); - auto tmp = analytical_solution_sphere_plan(alpha, 0., dt * (total_it + 1), radius, g, h); + auto tmp = analytical_solution_sphere_plane(alpha, 0., dt * (total_it + 1), radius, g, h); auto pos_analytical = tmp.first; auto q_analytical = quaternion(tmp.second); double error_pos = xt::linalg::norm(pos(1) - pos_analytical) / xt::linalg::norm(pos_analytical); double error_q = xt::linalg::norm(q(1) - q_analytical) / xt::linalg::norm(q_analytical); auto v = particles.v(); auto omega = particles.omega(); - tmp = analytical_solution_sphere_plan_velocity(alpha, 0., dt * (total_it + 1), radius, g, h); + tmp = analytical_solution_sphere_plane_velocity(alpha, 0., dt * (total_it + 1), radius, g, h); auto v_analytical = tmp.first; double error_v = xt::linalg::norm(v(1) - v_analytical) / xt::linalg::norm(v_analytical); @@ -307,11 +307,11 @@ namespace scopi REQUIRE(omega(1) == doctest::Approx(0.)); } - TEST_CASE_TEMPLATE_APPLY(sphere_plan, solver_dry_without_friction_t<2>); - TEST_CASE_TEMPLATE_APPLY(sphere_plan_force, solver_dry_without_friction_t<2, vap_fpd>); + TEST_CASE_TEMPLATE_APPLY(sphere_plane, solver_dry_without_friction_t<2>); + TEST_CASE_TEMPLATE_APPLY(sphere_plane_force, solver_dry_without_friction_t<2, vap_fpd>); TEST_CASE_TEMPLATE_APPLY(sphere_sphere_fixed, solver_dry_without_friction_t<2>); TEST_CASE_TEMPLATE_APPLY(sphere_sphere_fixed_force, solver_dry_without_friction_t<2, vap_fpd>); TEST_CASE_TEMPLATE_APPLY(sphere_sphere_moving, solver_dry_without_friction_t<2, vap_fpd>); - TEST_CASE_TEMPLATE_APPLY(sphere_inclined_plan, solver_dry_without_friction_t<2, vap_fpd>); + TEST_CASE_TEMPLATE_APPLY(sphere_inclined_plane, solver_dry_without_friction_t<2, vap_fpd>); } diff --git a/test/test_sphere.cpp b/test/test_sphere.cpp index eb49e704..44a49b3d 100644 --- a/test/test_sphere.cpp +++ b/test/test_sphere.cpp @@ -366,10 +366,10 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = 1; // total_it-1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_spheres_asymmetrical.json", tolerance)); } @@ -402,10 +402,10 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = 1; // total_it-1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_spheres_symmetrical.json", tolerance)); } @@ -459,11 +459,11 @@ namespace scopi } } - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = 1; // total_it-1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0019.json", "../test/references/2d_case_spheres.json", tolerance)); } diff --git a/test/test_superellipsoid.cpp b/test/test_superellipsoid.cpp index 948b655c..ba3cfd99 100644 --- a/test/test_superellipsoid.cpp +++ b/test/test_superellipsoid.cpp @@ -482,10 +482,10 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0199.json", "../test/references/two_ellipsoids_symmetrical.json", tolerance)); } @@ -522,10 +522,10 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0049.json", "../test/references/two_ellipsoids_spheres_symmetrical.json", tolerance)); } @@ -562,10 +562,10 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_ellipsoids_asymmetrical.json", tolerance)); } @@ -604,10 +604,10 @@ namespace scopi {-0.25, 0} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_ellipsoids_spheres_asymmetrical.json", tolerance)); } @@ -673,10 +673,10 @@ namespace scopi } } - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_ellipsoids_spheres_asymmetrical.json", tolerance)); } diff --git a/test/test_viscosity.cpp b/test/test_viscosity.cpp index 888e3618..9946206c 100644 --- a/test/test_viscosity.cpp +++ b/test/test_viscosity.cpp @@ -48,7 +48,7 @@ namespace scopi // set_params_test_uzawa(params); // } - TEST_CASE_TEMPLATE_DEFINE("sphere plan viscosity", SolverType, sphere_plan_viscosity) + TEST_CASE_TEMPLATE_DEFINE("sphere plane viscosity", SolverType, sphere_plane_viscosity) { constexpr std::size_t dim = 2; @@ -61,7 +61,7 @@ namespace scopi double dt = 0.05; scopi_container particles; - plan p( + plane p( { {0., 0.} }, @@ -77,20 +77,20 @@ namespace scopi {g * std::cos(alpha), -g * std::sin(alpha)} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); set_params_test(params.optim_params); params.solver_params.output_frequency = 1; - solver.run(150); + solver.run(dt, 150); particles.f()(1)(1) *= -1.; - solver.run(189, 150); + solver.run(dt, 189, 150); - CHECK(diffFile("./Results/scopi_objects_0188.json", "../test/references/sphere_plan_viscosity.json", tolerance)); + CHECK(diffFile("./Results/scopi_objects_0188.json", "../test/references/sphere_plane_viscosity.json", tolerance)); } - TEST_CASE_TEMPLATE_APPLY(sphere_plan_viscosity, solver_dry_without_friction_t<2>); + TEST_CASE_TEMPLATE_APPLY(sphere_plane_viscosity, solver_dry_without_friction_t<2>); - TEST_CASE_TEMPLATE_DEFINE("sphere plan viscosity friction vertical", SolverType, sphere_plan_viscosity_friction_vertical) + TEST_CASE_TEMPLATE_DEFINE("sphere plane viscosity friction vertical", SolverType, sphere_plane_viscosity_friction_vertical) { constexpr std::size_t dim = 2; @@ -103,7 +103,7 @@ namespace scopi std::size_t total_it = 100; scopi_container particles; - plan p( + plane p( { {0., 0.} }, @@ -119,20 +119,20 @@ namespace scopi {0, -g} })); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); // params.problem_params.mu = 0.1; params.solver_params.output_frequency = 1; - solver.run(total_it); + solver.run(dt, total_it); particles.f()(1)(1) *= -1.; - solver.run(2 * total_it, total_it); + solver.run(dt, 2 * total_it, total_it); - CHECK(diffFile("./Results/scopi_objects_0199.json", "../test/references/sphere_plan_viscosity_friction_vertical.json", tolerance)); + CHECK(diffFile("./Results/scopi_objects_0199.json", "../test/references/sphere_plane_viscosity_friction_vertical.json", tolerance)); } - TEST_CASE_TEMPLATE_APPLY(sphere_plan_viscosity_friction_vertical, solver_dry_with_friction_t<2, vap_fpd>); + TEST_CASE_TEMPLATE_APPLY(sphere_plane_viscosity_friction_vertical, solver_dry_with_friction_t<2, vap_fpd>); - TEST_CASE_TEMPLATE_DEFINE("sphere plan viscosity friction", SolverType, sphere_plan_viscosity_friction) + TEST_CASE_TEMPLATE_DEFINE("sphere plane viscosity friction", SolverType, sphere_plane_viscosity_friction) { constexpr std::size_t dim = 2; @@ -145,7 +145,7 @@ namespace scopi std::size_t total_it = 100; scopi_container particles; - plan p( + plane p( { {0., 0.} }, @@ -161,16 +161,16 @@ namespace scopi {g, -g} })); - SolverType solver(particles, dt); + SolverType solver(particles); // auto params = solver.get_params(); // params.problem_params.mu = 0.1; - solver.run(total_it); + solver.run(dt, total_it); particles.f()(1)(1) *= -1.; - solver.run(2 * total_it, total_it); + solver.run(dt, 2 * total_it, total_it); - CHECK(diffFile("./Results/scopi_objects_0199.json", "../test/references/sphere_plan_viscosity_friction.json", tolerance)); + CHECK(diffFile("./Results/scopi_objects_0199.json", "../test/references/sphere_plane_viscosity_friction.json", tolerance)); } - TEST_CASE_TEMPLATE_APPLY(sphere_plan_viscosity_friction, solver_dry_with_friction_t<2, vap_fpd>); + TEST_CASE_TEMPLATE_APPLY(sphere_plane_viscosity_friction, solver_dry_with_friction_t<2, vap_fpd>); } diff --git a/test/test_worm.cpp b/test/test_worm.cpp index 4ef2657c..4c4bd206 100644 --- a/test/test_worm.cpp +++ b/test/test_worm.cpp @@ -512,12 +512,12 @@ namespace scopi particles.push_back(w1, prop.desired_velocity({-1., 0.})); particles.push_back(w2, prop.desired_velocity({1., 0.})); - SolverType solver(particles, dt); + SolverType solver(particles); auto params = solver.get_params(); set_params_test(params.optim_params); params.contact_method_params.dmax = 1.; // params.solver_params.output_frequency = total_it - 1; - solver.run(total_it); + solver.run(dt, total_it); CHECK(diffFile("./Results/scopi_objects_0999.json", "../test/references/two_worms.json", tolerance)); } diff --git a/test/utils.cpp b/test/utils.cpp index 37349867..6c5d8f7f 100644 --- a/test/utils.cpp +++ b/test/utils.cpp @@ -57,7 +57,7 @@ namespace scopi return false; } - std::pair, double> analytical_solution_sphere_plan(double alpha, double mu, double t, double r, double g, double y0) + std::pair, double> analytical_solution_sphere_plane(double alpha, double mu, double t, double r, double g, double y0) { const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); type::position_t<2> x; @@ -90,7 +90,7 @@ namespace scopi } std::pair, double> - analytical_solution_sphere_plan_velocity(double alpha, double mu, double t, double r, double g, double y0) + analytical_solution_sphere_plane_velocity(double alpha, double mu, double t, double r, double g, double y0) { const double t_impact = std::sqrt(2 * (y0 - r) / (g * std::cos(alpha))); type::position_t<2> x; diff --git a/test/utils.hpp b/test/utils.hpp index 2819d926..eae84a4b 100644 --- a/test/utils.hpp +++ b/test/utils.hpp @@ -11,7 +11,7 @@ namespace scopi static constexpr double PI = xt::numeric_constants::PI; static constexpr double tolerance = 1e-6; bool diffFile(const std::string& filenameRef, const std::string& filenameResult, double tol = 1e-12); - std::pair, double> analytical_solution_sphere_plan(double alpha, double mu, double t, double r, double g, double y0); + std::pair, double> analytical_solution_sphere_plane(double alpha, double mu, double t, double r, double g, double y0); std::pair, double> - analytical_solution_sphere_plan_velocity(double alpha, double mu, double t, double r, double g, double y0); + analytical_solution_sphere_plane_velocity(double alpha, double mu, double t, double r, double g, double y0); } diff --git a/tools/visualization/js/src/index.js b/tools/visualization/js/src/index.js index 8cff6190..518712fc 100644 --- a/tools/visualization/js/src/index.js +++ b/tools/visualization/js/src/index.js @@ -16,6 +16,7 @@ var options = { "pause": true, "contacts" : true, "radiuses" : true, + "show_plan" : false, "current_frame": 0, }; @@ -176,7 +177,7 @@ function drawObjects() { nbSpheres = 0; objects.forEach((obj, index) => { - if (obj.type === "plan") { + if (options.show_plan && obj.type === "plane") { const normal = new THREE.Vector3(); const position = new THREE.Vector3(); // TODO this is 3D only @@ -205,7 +206,7 @@ function drawObjects() { nbSpheres += 1; }); } - else if (obj.type === "sphere") { + else if (obj.type === "sphere" || obj.type === "superellipsoid") { sphereObject(obj, matrix, rot); mesh.setMatrixAt(nbSpheres, matrix); nbSpheres += 1; @@ -403,6 +404,10 @@ function init() { function (value) { drawObjects(); }); + view.add(options, 'show_plan').name('Show plans').listen().onChange( + function (value) { + drawObjects(); + }); view.open(); var obj = { diff --git a/tools/visualization/python/test_visu.py b/tools/visualization/python/test_visu.py index 07bc11f6..c78c0295 100644 --- a/tools/visualization/python/test_visu.py +++ b/tools/visualization/python/test_visu.py @@ -106,7 +106,7 @@ geom["e"] = obj["squareness"][0]*np.ones((np.array(geom.points).shape[0],)) geom["n"] = obj["squareness"][1]*np.ones((np.array(geom.points).shape[0],)) geom["e+n"] = obj["squareness"][0]*np.ones((np.array(geom.points).shape[0],))+obj["squareness"][1]*np.ones((np.array(geom.points).shape[0],)) - elif (obj["type"] == "plan"): + elif (obj["type"] == "plane"): if (len(obj["position"])==2): # 2D geom = pv.Plane( center=(obj["position"][0],obj["position"][1],0), @@ -187,7 +187,7 @@ geom["n"] = obj["squareness"][1]*np.ones((np.array(geom.points).shape[0],)) geom["e+n"] = obj["squareness"][0]*np.ones((np.array(geom.points).shape[0],))+obj["squareness"][1]*np.ones((np.array(geom.points).shape[0],)) - elif (obj["type"] == "plan"): + elif (obj["type"] == "plane"): if (len(obj["position"])==2): # 2D geom = pv.Plane( center=(obj["position"][0],obj["position"][1],0),