diff --git a/src/DEM/API.h b/src/DEM/API.h index c1ff1665..c111d0a6 100644 --- a/src/DEM/API.h +++ b/src/DEM/API.h @@ -1752,11 +1752,11 @@ class DEMSolver { // Some float3 quantity that is representitive of a component's initial orientation (such as plane normal, and its // meaning can vary among different types) std::vector m_anal_comp_rot; - // Some float quantity that is representitive of a component's size (e.g. for a cylinder, top radius) + // Some float quantity that is representitive of a component's size (e.g. cylinder radius or cone slope) std::vector m_anal_size_1; - // Some float quantity that is representitive of a component's size (e.g. for a cylinder, bottom radius) + // Some float quantity that is representitive of a component's size (e.g. cone lower axial bound) std::vector m_anal_size_2; - // Some float quantity that is representitive of a component's size (e.g. for a cylinder, its length) + // Some float quantity that is representitive of a component's size (e.g. cone upper axial bound) std::vector m_anal_size_3; // Component object types std::vector m_anal_types; diff --git a/src/DEM/APIPrivate.cpp b/src/DEM/APIPrivate.cpp index 08250829..e3a24abf 100644 --- a/src/DEM/APIPrivate.cpp +++ b/src/DEM/APIPrivate.cpp @@ -683,6 +683,16 @@ void DEMSolver::preprocessAnalyticalObjs() { addAnalCompTemplate(ANAL_OBJ_TYPE_CYL_INF, comp_mat.at(i), thisExtObj, param.cyl.center, param.cyl.dir, param.cyl.radius, 0, 0, param.cyl.normal); break; + case OBJ_COMPONENT::CONE_INF: + addAnalCompTemplate(ANAL_OBJ_TYPE_CONE_INF, comp_mat.at(i), thisExtObj, param.cone.cone_tip, + param.cone.dir, param.cone.slope, param.cone.hmin, param.cone.hmax, + param.cone.normal); + break; + case OBJ_COMPONENT::CONE: + addAnalCompTemplate(ANAL_OBJ_TYPE_CONE, comp_mat.at(i), thisExtObj, param.cone.cone_tip, + param.cone.dir, param.cone.slope, param.cone.hmin, param.cone.hmax, + param.cone.normal); + break; default: DEME_ERROR("There is at least one analytical boundary that has a type not supported."); } diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index 98473bf0..23d589e6 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -33,9 +33,11 @@ struct DEMSphereParams_t { objNormal_t normal; }; -/// Cone (pos is tip pos) +/// Cone side. `cone_tip' is the apex location, `dir' points from the apex toward the base, and `slope' is radius +/// growth per unit axial distance. struct DEMConeParams_t { float3 cone_tip; + float3 dir; float slope; float hmax; float hmin; @@ -64,6 +66,29 @@ struct DEMCylinderParams_t { objNormal_t normal; }; +inline void assertConeInputs(const float3 axis, + const float slope, + const float hmin, + const float hmax, + const std::string& func_name) { + if (!std::isfinite(axis.x) || !std::isfinite(axis.y) || !std::isfinite(axis.z) || + length(axis) <= DEME_TINY_FLOAT) { + std::stringstream out; + out << func_name << "'s axis argument must be a finite, non-zero direction.\n"; + throw std::runtime_error(out.str()); + } + if (!std::isfinite(slope) || slope <= 0.0f) { + std::stringstream out; + out << func_name << "'s slope argument must be finite and positive.\n"; + throw std::runtime_error(out.str()); + } + if (!std::isfinite(hmin) || !std::isfinite(hmax) || hmin < 0.0f || hmax <= hmin) { + std::stringstream out; + out << func_name << "'s axial bounds must be finite and satisfy 0 <= hmin < hmax.\n"; + throw std::runtime_error(out.str()); + } +} + /// API-(Host-)side struct that holds cached user-input external objects class DEMExternObj : public DEMInitializer { public: @@ -99,6 +124,7 @@ class DEMExternObj : public DEMInitializer { DEMPlateParams_t plate; DEMPlaneParams_t plane; DEMCylinderParams_t cyl; + DEMConeParams_t cone; }; std::vector entity_params; @@ -224,6 +250,74 @@ class DEMExternObj : public DEMInitializer { assertThreeElements(axis, "AddCylinder", "axis"); AddCylinder(make_float3(pos[0], pos[1], pos[2]), make_float3(axis[0], axis[1], axis[2]), rad, material, normal); } + + /// Add a single-nappe cone side that starts at `tip' and extends indefinitely along `axis'. + /// + /// The cone is an analytical shell, like analytical cylinders. `slope' is radius per unit axial distance from the + /// tip. Use ENTITY_NORMAL_INWARD for particles inside the cone and ENTITY_NORMAL_OUTWARD for particles outside a + /// solid cone. + void AddCone(const float3 tip, + const float3 axis, + const float slope, + const std::shared_ptr& material, + const objNormal_t normal = ENTITY_NORMAL_INWARD) { + assertConeInputs(axis, slope, 0.0f, DEME_HUGE_FLOAT, "AddCone"); + types.push_back(OBJ_COMPONENT::CONE_INF); + materials.push_back(material); + DEMAnalEntParams params; + params.cone.cone_tip = tip; + params.cone.dir = normalize(axis); + params.cone.slope = slope; + params.cone.hmin = 0.0f; + params.cone.hmax = DEME_HUGE_FLOAT; + params.cone.normal = normal; + entity_params.push_back(params); + } + void AddCone(const std::vector& tip, + const std::vector& axis, + const float slope, + const std::shared_ptr& material, + const objNormal_t normal = ENTITY_NORMAL_INWARD) { + assertThreeElements(tip, "AddCone", "tip"); + assertThreeElements(axis, "AddCone", "axis"); + AddCone(make_float3(tip[0], tip[1], tip[2]), make_float3(axis[0], axis[1], axis[2]), slope, material, normal); + } + + /// Add a height-bounded cone side. The surface is clipped to hmin <= dot(P - tip, axis) <= hmax. + /// + /// This creates the exact conical side/rim shell only; add analytical planes separately when cap contacts are + /// needed. Setting hmin > 0 gives a conical frustum side without introducing mesh facets. + void AddConeSegment(const float3 tip, + const float3 axis, + const float slope, + const float hmin, + const float hmax, + const std::shared_ptr& material, + const objNormal_t normal = ENTITY_NORMAL_INWARD) { + assertConeInputs(axis, slope, hmin, hmax, "AddConeSegment"); + types.push_back(OBJ_COMPONENT::CONE); + materials.push_back(material); + DEMAnalEntParams params; + params.cone.cone_tip = tip; + params.cone.dir = normalize(axis); + params.cone.slope = slope; + params.cone.hmin = hmin; + params.cone.hmax = hmax; + params.cone.normal = normal; + entity_params.push_back(params); + } + void AddConeSegment(const std::vector& tip, + const std::vector& axis, + const float slope, + const float hmin, + const float hmax, + const std::shared_ptr& material, + const objNormal_t normal = ENTITY_NORMAL_INWARD) { + assertThreeElements(tip, "AddConeSegment", "tip"); + assertThreeElements(axis, "AddConeSegment", "axis"); + AddConeSegment(make_float3(tip[0], tip[1], tip[2]), make_float3(axis[0], axis[1], axis[2]), slope, hmin, hmax, + material, normal); + } }; // DEM mesh object diff --git a/src/DEM/Defines.h b/src/DEM/Defines.h index c2a8de8e..7220c5ae 100644 --- a/src/DEM/Defines.h +++ b/src/DEM/Defines.h @@ -68,6 +68,8 @@ constexpr clumpComponentOffset_t NUM_ACTIVE_TEMPLATE_LOADING_THREADS = const objType_t ANAL_OBJ_TYPE_PLANE = 0; const objType_t ANAL_OBJ_TYPE_PLATE = 1; const objType_t ANAL_OBJ_TYPE_CYL_INF = 2; +const objType_t ANAL_OBJ_TYPE_CONE_INF = 3; +const objType_t ANAL_OBJ_TYPE_CONE = 4; const objNormal_t ENTITY_NORMAL_INWARD = 0; const objNormal_t ENTITY_NORMAL_OUTWARD = 1; diff --git a/src/demo/CMakeLists.txt b/src/demo/CMakeLists.txt index 9aca34a9..f529b0a8 100644 --- a/src/demo/CMakeLists.txt +++ b/src/demo/CMakeLists.txt @@ -22,6 +22,8 @@ SET(DEMOS DEMdemo_TestPack DEMdemo_RotatingDrum DEMdemo_Centrifuge + DEMdemo_AnalyticalCone + DEMdemo_AnalyticalConePour DEMdemo_GameOfLife DEMdemo_BallDrop DEMdemo_BallDrop2D @@ -90,4 +92,3 @@ FOREACH(PROGRAM ${DEMOS}) # install(TARGETS ${PROGRAM} DESTINATION ${DEME_INSTALL_DEMO}) ENDFOREACH(PROGRAM) - diff --git a/src/demo/DEMdemo_AnalyticalCone.cpp b/src/demo/DEMdemo_AnalyticalCone.cpp new file mode 100644 index 00000000..d68b0648 --- /dev/null +++ b/src/demo/DEMdemo_AnalyticalCone.cpp @@ -0,0 +1,210 @@ +// Copyright (c) 2021, SBEL GPU Development Team +// Copyright (c) 2021, University of Wisconsin - Madison +// +// SPDX-License-Identifier: BSD-3-Clause + +// ============================================================================= +// This small validation demo exercises analytical cone contacts without using +// triangle meshes. It keeps all contacts in one DEM solver initialization: +// - inward cone side contact, +// - outward cone side contact, +// - apex contact for the single-nappe cone tip, +// - a no-contact control. +// ============================================================================= + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace deme; + +namespace { + +struct ConeValidationCase { + std::string name; + std::string title; + float3 sphere_position; + float3 cone_tip; + float3 cone_axis; + float slope; + objNormal_t normal; + unsigned int sphere_family; + unsigned int cone_family; + int expected_force_x_sign; + int expected_force_z_sign; + bool expected_active_contact; +}; + +float length3(const float3& v) { + return std::sqrt(v.x * v.x + v.y * v.y + v.z * v.z); +} + +void require_condition(bool condition, const std::string& message) { + if (!condition) { + throw std::runtime_error(message); + } +} + +void require_force_direction(const std::map& forces, + bodyID_t owner, + const std::string& case_name, + int x_sign, + int z_sign) { + auto search = forces.find(owner); + require_condition(search != forces.end(), case_name + " did not produce an active contact."); + + const float3 force = search->second; + require_condition(length3(force) > 0.0f, case_name + " contact force is zero."); + if (x_sign < 0) { + require_condition(force.x < 0.0f, case_name + " force.x should be negative."); + } else if (x_sign > 0) { + require_condition(force.x > 0.0f, case_name + " force.x should be positive."); + } else { + require_condition(std::fabs(force.x) < 1.0e-3f * length3(force), case_name + " force.x should be near zero."); + } + + if (z_sign < 0) { + require_condition(force.z < 0.0f, case_name + " force.z should be negative."); + } else if (z_sign > 0) { + require_condition(force.z > 0.0f, case_name + " force.z should be positive."); + } else { + require_condition(std::fabs(force.z) < 1.0e-3f * length3(force), case_name + " force.z should be near zero."); + } +} + +const char* normal_mode_name(objNormal_t normal) { + return (normal == ENTITY_NORMAL_OUTWARD) ? "outward" : "inward"; +} + +void write_validation_case_metadata(const std::filesystem::path& out_file, + const std::vector& cases, + float sphere_radius) { + std::ofstream metadata(out_file); + require_condition(static_cast(metadata), "Unable to open validation case metadata file."); + + metadata << "case_id,name,title,sphere_id,sphere_family,cone_family,sphere_radius," + << "sphere_initial_x,sphere_initial_y,sphere_initial_z," + << "tip_x,tip_y,tip_z,axis_x,axis_y,axis_z,slope,normal_mode," + << "expected_active_contact,expected_force_x_sign,expected_force_z_sign,display_height\n"; + for (size_t i = 0; i < cases.size(); i++) { + const auto& this_case = cases.at(i); + // The display height only controls the ParaView inspection shell. The + // DEM object itself remains the semi-infinite analytical AddCone used + // below, so this metadata does not introduce a mesh into the test. + constexpr float display_height = 2.25f; + metadata << i << "," << this_case.name << "," << this_case.title << "," << i << "," + << this_case.sphere_family << "," << this_case.cone_family << "," << sphere_radius << "," + << this_case.sphere_position.x << "," << this_case.sphere_position.y << "," + << this_case.sphere_position.z << "," << this_case.cone_tip.x << "," << this_case.cone_tip.y << "," + << this_case.cone_tip.z << "," << this_case.cone_axis.x << "," << this_case.cone_axis.y << "," + << this_case.cone_axis.z << "," << this_case.slope << "," << normal_mode_name(this_case.normal) << "," + << (this_case.expected_active_contact ? 1 : 0) << "," << this_case.expected_force_x_sign << "," + << this_case.expected_force_z_sign << "," << display_height << "\n"; + } +} + +} // namespace + +int main(int argc, char** argv) { + const std::filesystem::path out_dir = + (argc > 1) ? std::filesystem::path(argv[1]) : std::filesystem::path("/tmp/deme_cone_validation"); + std::filesystem::create_directories(out_dir); + + DEMSolver DEMSim(1); + DEMSim.SetVerbosity(INFO); + DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV); + DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV); + DEMSim.SetContactOutputContent({"OWNER", "GEO_ID", "FORCE", "POINT", "NORMAL"}); + DEMSim.UseFrictionlessHertzianModel(); + + auto mat_sphere = DEMSim.LoadMaterial({{"E", 1.0e7f}, {"nu", 0.3f}, {"CoR", 0.1f}, {"mu", 0.0f}, {"Crr", 0.0f}}); + auto mat_cone = DEMSim.LoadMaterial({{"E", 1.0e7f}, {"nu", 0.3f}, {"CoR", 0.1f}, {"mu", 0.0f}, {"Crr", 0.0f}}); + + constexpr float sphere_radius = 0.2f; + auto sphere_type = DEMSim.LoadSphereType(1.0f, sphere_radius, mat_sphere); + + const std::vector cases = { + {"inward_side", "Inward side contact", make_float3(0.9f, 0.0f, 1.0f), make_float3(0.0f, 0.0f, 0.0f), + make_float3(0.0f, 0.0f, 1.0f), 1.0f, ENTITY_NORMAL_INWARD, 10, 20, -1, 1, true}, + {"outward_side", "Outward side contact", make_float3(5.1f, 0.0f, 1.0f), make_float3(4.0f, 0.0f, 0.0f), + make_float3(0.0f, 0.0f, 1.0f), 1.0f, ENTITY_NORMAL_OUTWARD, 11, 21, 1, -1, true}, + {"apex", "Apex contact", make_float3(8.0f, 0.0f, -0.05f), make_float3(8.0f, 0.0f, 0.0f), + make_float3(0.0f, 0.0f, 1.0f), 1.0f, ENTITY_NORMAL_INWARD, 12, 22, 0, -1, true}, + {"no_contact", "No-contact control", make_float3(12.5f, 0.0f, 1.0f), make_float3(12.0f, 0.0f, 0.0f), + make_float3(0.0f, 0.0f, 1.0f), 1.0f, ENTITY_NORMAL_INWARD, 13, 23, 0, 0, false} + }; + + std::vector sphere_positions; + std::vector sphere_families; + for (const auto& this_case : cases) { + sphere_positions.push_back(this_case.sphere_position); + sphere_families.push_back(this_case.sphere_family); + } + auto spheres = DEMSim.AddClumps(sphere_type, sphere_positions); + spheres->SetFamilies(sphere_families); + + for (const auto& this_case : cases) { + auto cone = DEMSim.AddExternalObject(); + cone->AddCone(this_case.cone_tip, this_case.cone_axis, this_case.slope, mat_cone, this_case.normal); + cone->SetFamily(this_case.cone_family); + } + + for (const auto& sphere_case : cases) { + for (const auto& cone_case : cases) { + if (sphere_case.cone_family != cone_case.cone_family) { + DEMSim.DisableContactBetweenFamilies(sphere_case.sphere_family, cone_case.cone_family); + } + } + } + + for (const auto& this_case : cases) { + DEMSim.SetFamilyFixed(this_case.cone_family); + } + DEMSim.InstructBoxDomainDimension({-1.5f, 14.0f}, {-1.5f, 1.5f}, {-1.0f, 2.5f}); + DEMSim.InstructBoxDomainBoundingBC("none", mat_cone); + DEMSim.SetInitTimeStep(1.0e-6); + DEMSim.SetCDUpdateFreq(1); + DEMSim.SetGravitationalAcceleration(make_float3(0.0f)); + DEMSim.SetMaxVelocity(1.0f); + + DEMSim.Initialize(); + DEMSim.DoDynamicsThenSync(1.0e-6); + + auto potential_contacts = DEMSim.GetContactDetailedInfo(-1.0f); + DEMSim.WriteSphereFile(out_dir / "spheres.csv"); + DEMSim.WriteContactFile((out_dir / "contacts_potential.csv").string(), -1.0f); + DEMSim.WriteContactFile((out_dir / "contacts_active.csv").string(), DEME_TINY_FLOAT); + write_validation_case_metadata(out_dir / "validation_cases.csv", cases, sphere_radius); + + std::cout << "Potential analytical contact entries: " << potential_contacts->Size() << std::endl; + + auto contacts = DEMSim.GetContactDetailedInfo(DEME_TINY_FLOAT); + auto& owners = contacts->GetAOwner(); + auto& forces = contacts->GetForce(); + require_condition(owners.size() == forces.size(), "Contact owner/force arrays have inconsistent lengths."); + std::cout << "Active analytical contact entries: " << forces.size() << std::endl; + require_condition(forces.size() == 3, "Expected exactly three active analytical cone contacts."); + + std::map force_by_owner; + for (size_t i = 0; i < forces.size(); i++) { + force_by_owner[owners[i]] = forces[i]; + std::cout << "owner=" << owners[i] << " force=(" << forces[i].x << ", " << forces[i].y << ", " << forces[i].z + << ") normal=(" << contacts->GetNormal()[i].x << ", " << contacts->GetNormal()[i].y << ", " + << contacts->GetNormal()[i].z << ") type=" << contacts->GetContactType()[i] << std::endl; + } + + require_force_direction(force_by_owner, 0, "inward cone side", -1, 1); + require_force_direction(force_by_owner, 1, "outward cone side", 1, -1); + require_force_direction(force_by_owner, 2, "cone apex", 0, -1); + require_condition(force_by_owner.find(3) == force_by_owner.end(), "no-contact control produced an active contact."); + + std::cout << "Analytical cone GPU validation passed. Output directory: " << out_dir << std::endl; + return 0; +} diff --git a/src/demo/DEMdemo_AnalyticalConePour.cpp b/src/demo/DEMdemo_AnalyticalConePour.cpp new file mode 100644 index 00000000..db594e62 --- /dev/null +++ b/src/demo/DEMdemo_AnalyticalConePour.cpp @@ -0,0 +1,249 @@ +// Copyright (c) 2021, SBEL GPU Development Team +// Copyright (c) 2021, University of Wisconsin - Madison +// +// SPDX-License-Identifier: BSD-3-Clause + +// ============================================================================= +// This demo pours spheres into three analytical cone/cylinder containers: +// 1. a cone with its corner down, used as a hollow conical wall, +// 2. a cone with its corner up, used as a solid conical obstacle, +// 3. a cone whose analytical tip is below the floor, so the active surface is +// a conical frustum. +// +// The cone and cylinder contacts are analytical DEM-Engine objects. The output +// CSV files can be converted into a ParaView time series for visual inspection. +// ============================================================================= + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace deme; + +namespace { + +struct PourCase { + std::string name; + std::string title; + float x_center; + float3 cone_tip; + float3 cone_axis; + float slope; + float hmin; + float hmax; + objNormal_t cone_normal; + float cylinder_radius; + float bottom_z; + float top_z; + unsigned int sphere_family; + unsigned int boundary_family; +}; + +struct SphereSeed { + size_t sphere_id; + unsigned int case_id; + float radius; + float3 position; + float3 velocity; +}; + +void require_condition(bool condition, const std::string& message) { + if (!condition) { + throw std::runtime_error(message); + } +} + +const char* normal_mode_name(objNormal_t normal) { + return (normal == ENTITY_NORMAL_OUTWARD) ? "outward" : "inward"; +} + +std::string frame_file_name(const std::string& stem, int frame) { + std::ostringstream name; + name << stem << "_" << std::setw(4) << std::setfill('0') << frame << ".csv"; + return name.str(); +} + +void write_pour_case_metadata(const std::filesystem::path& out_file, const std::vector& cases) { + std::ofstream metadata(out_file); + require_condition(static_cast(metadata), "Unable to open pour case metadata file."); + + metadata << "case_id,name,title,x_center,tip_x,tip_y,tip_z,axis_x,axis_y,axis_z,slope,hmin,hmax," + << "normal_mode,cylinder_radius,bottom_z,top_z,sphere_family,boundary_family\n"; + for (size_t i = 0; i < cases.size(); i++) { + const auto& this_case = cases.at(i); + metadata << i << "," << this_case.name << "," << this_case.title << "," << this_case.x_center << "," + << this_case.cone_tip.x << "," << this_case.cone_tip.y << "," << this_case.cone_tip.z << "," + << this_case.cone_axis.x << "," << this_case.cone_axis.y << "," << this_case.cone_axis.z << "," + << this_case.slope << "," << this_case.hmin << "," << this_case.hmax << "," + << normal_mode_name(this_case.cone_normal) << "," << this_case.cylinder_radius << "," + << this_case.bottom_z << "," << this_case.top_z << "," << this_case.sphere_family << "," + << this_case.boundary_family << "\n"; + } +} + +void write_pour_sphere_metadata(const std::filesystem::path& out_file, const std::vector& spheres) { + std::ofstream metadata(out_file); + require_condition(static_cast(metadata), "Unable to open pour sphere metadata file."); + + metadata << "sphere_id,case_id,radius,initial_x,initial_y,initial_z,initial_vx,initial_vy,initial_vz\n"; + for (const auto& sphere : spheres) { + metadata << sphere.sphere_id << "," << sphere.case_id << "," << sphere.radius << "," << sphere.position.x + << "," << sphere.position.y << "," << sphere.position.z << "," << sphere.velocity.x << "," + << sphere.velocity.y << "," << sphere.velocity.z << "\n"; + } +} + +void write_dynamic_frame_outputs(DEMSolver& DEMSim, const std::filesystem::path& dynamic_dir, int frame) { + DEMSim.WriteSphereFile(dynamic_dir / frame_file_name("spheres", frame)); + DEMSim.WriteContactFile((dynamic_dir / frame_file_name("contacts_active", frame)).string(), DEME_TINY_FLOAT); +} + +std::vector make_sphere_seeds(const std::vector& cases, float radius) { + std::vector seeds; + constexpr float spacing = 0.118f; + constexpr float downward_speed = -0.95f; + + for (size_t case_id = 0; case_id < cases.size(); case_id++) { + const auto& this_case = cases.at(case_id); + for (int layer = 0; layer < 3; layer++) { + for (int ix = -2; ix <= 2; ix++) { + for (int iy = -2; iy <= 2; iy++) { + const float dx = spacing * ix + 0.018f * ((layer + iy + 8) % 3 - 1); + const float dy = spacing * iy + 0.014f * ((layer + ix + 8) % 3 - 1); + if (std::sqrt(dx * dx + dy * dy) > this_case.cylinder_radius - 3.2f * radius) { + continue; + } + const float z = + this_case.top_z + 0.16f + 0.12f * layer + 0.012f * static_cast((ix + iy + 8) % 3); + seeds.push_back({seeds.size(), static_cast(case_id), radius, + make_float3(this_case.x_center + dx, dy, z), + make_float3(0.0f, 0.0f, downward_speed)}); + } + } + } + } + return seeds; +} + +} // namespace + +int main(int argc, char** argv) { + const std::filesystem::path out_dir = + (argc > 1) ? std::filesystem::path(argv[1]) : std::filesystem::path("/tmp/deme_cone_pour_validation"); + const std::filesystem::path dynamic_dir = out_dir / "dynamic"; + std::filesystem::create_directories(dynamic_dir); + + DEMSolver DEMSim(1); + DEMSim.SetVerbosity(INFO); + DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV); + DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV); + DEMSim.SetContactOutputContent({"OWNER", "GEO_ID", "FORCE", "POINT", "NORMAL"}); + DEMSim.UseFrictionalHertzianModel(); + + auto mat_sphere = + DEMSim.LoadMaterial({{"E", 2.0e5f}, {"nu", 0.3f}, {"CoR", 0.2f}, {"mu", 0.45f}, {"Crr", 0.01f}}); + auto mat_wall = + DEMSim.LoadMaterial({{"E", 2.0e5f}, {"nu", 0.3f}, {"CoR", 0.2f}, {"mu", 0.45f}, {"Crr", 0.01f}}); + + constexpr float cylinder_radius = 0.72f; + constexpr float bottom_z = 0.0f; + constexpr float cone_height = 1.25f; + constexpr float top_z = bottom_z + cone_height; + constexpr float frustum_tip_below_floor = 1.0f; + const float frustum_slope = cylinder_radius / (frustum_tip_below_floor + cone_height); + + const std::vector cases = { + {"corner_down", "Cone corner down", -2.25f, make_float3(-2.25f, 0.0f, bottom_z), + make_float3(0.0f, 0.0f, 1.0f), cylinder_radius / cone_height, 0.0f, cone_height, ENTITY_NORMAL_INWARD, + cylinder_radius, bottom_z, top_z, 10, 110}, + {"corner_up", "Cone corner up", 0.0f, make_float3(0.0f, 0.0f, top_z), make_float3(0.0f, 0.0f, -1.0f), + cylinder_radius / cone_height, 0.0f, cone_height, ENTITY_NORMAL_OUTWARD, cylinder_radius, bottom_z, top_z, 11, + 111}, + {"frustum_from_below", "Cone tip below floor visible frustum", 2.25f, + make_float3(2.25f, 0.0f, bottom_z - frustum_tip_below_floor), make_float3(0.0f, 0.0f, 1.0f), + frustum_slope, frustum_tip_below_floor, frustum_tip_below_floor + cone_height, ENTITY_NORMAL_INWARD, + cylinder_radius, bottom_z, top_z, 12, 112}, + }; + + for (const auto& this_case : cases) { + auto container = DEMSim.AddExternalObject(); + container->AddCylinder(make_float3(this_case.x_center, 0.0f, bottom_z), make_float3(0.0f, 0.0f, 1.0f), + this_case.cylinder_radius, mat_wall, ENTITY_NORMAL_INWARD); + container->AddPlane(make_float3(this_case.x_center, 0.0f, this_case.bottom_z), make_float3(0.0f, 0.0f, 1.0f), + mat_wall); + container->AddConeSegment(this_case.cone_tip, this_case.cone_axis, this_case.slope, this_case.hmin, + this_case.hmax, mat_wall, this_case.cone_normal); + container->SetFamily(this_case.boundary_family); + DEMSim.SetFamilyFixed(this_case.boundary_family); + } + + constexpr float sphere_radius = 0.045f; + constexpr float sphere_density = 1800.0f; + constexpr float sphere_volume = 4.0f / 3.0f * 3.14159265358979323846f * sphere_radius * sphere_radius * sphere_radius; + const float sphere_mass = sphere_density * sphere_volume; + auto sphere_type = DEMSim.LoadSphereType(sphere_mass, sphere_radius, mat_sphere); + const std::vector seeds = make_sphere_seeds(cases, sphere_radius); + + std::vector sphere_positions; + std::vector sphere_velocities; + std::vector sphere_families; + for (const auto& seed : seeds) { + sphere_positions.push_back(seed.position); + sphere_velocities.push_back(seed.velocity); + sphere_families.push_back(cases.at(seed.case_id).sphere_family); + } + auto spheres = DEMSim.AddClumps(sphere_type, sphere_positions); + spheres->SetFamilies(sphere_families); + spheres->SetVel(sphere_velocities); + + for (const auto& sphere_case : cases) { + for (const auto& boundary_case : cases) { + if (sphere_case.boundary_family != boundary_case.boundary_family) { + DEMSim.DisableContactBetweenFamilies(sphere_case.sphere_family, boundary_case.boundary_family); + } + if (sphere_case.sphere_family != boundary_case.sphere_family) { + DEMSim.DisableContactBetweenFamilies(sphere_case.sphere_family, boundary_case.sphere_family); + } + } + } + + DEMSim.InstructBoxDomainDimension({-3.35f, 3.35f}, {-0.95f, 0.95f}, {-0.20f, 2.20f}); + DEMSim.InstructBoxDomainBoundingBC("none", mat_wall); + DEMSim.SetInitTimeStep(2.0e-5); + DEMSim.SetCDUpdateFreq(10); + DEMSim.SetGravitationalAcceleration(make_float3(0.0f, 0.0f, -9.81f)); + DEMSim.SetMaxVelocity(80.0f); + DEMSim.SetErrorOutVelocity(100.0f); + + write_pour_case_metadata(out_dir / "pour_cases.csv", cases); + write_pour_sphere_metadata(out_dir / "pour_spheres.csv", seeds); + + DEMSim.Initialize(); + + std::ofstream frame_times(dynamic_dir / "frame_times.csv"); + require_condition(static_cast(frame_times), "Unable to open dynamic frame time file."); + frame_times << "frame,time\n"; + + constexpr int dynamic_frames = 250; + constexpr double dynamic_frame_time = 2.0e-2; + for (int frame = 0; frame < dynamic_frames; frame++) { + frame_times << frame << "," << DEMSim.GetSimTime() << "\n"; + write_dynamic_frame_outputs(DEMSim, dynamic_dir, frame); + DEMSim.DoDynamicsThenSync(dynamic_frame_time); + if (frame % 30 == 0) { + std::cout << "Pour frame " << frame << " / " << dynamic_frames << std::endl; + } + } + + std::cout << "Analytical cone pouring validation wrote " << dynamic_frames << " frames for " << seeds.size() + << " spheres in " << out_dir << std::endl; + return 0; +} diff --git a/src/demo/DEMdemo_ConePenetration.cpp b/src/demo/DEMdemo_ConePenetration.cpp index 4931b603..2dcc838b 100644 --- a/src/demo/DEMdemo_ConePenetration.cpp +++ b/src/demo/DEMdemo_ConePenetration.cpp @@ -78,26 +78,26 @@ int main() { DEMSim.AddClumps(my_template, input_xyz); std::cout << "Total num of particles: " << input_xyz.size() << std::endl; - // Load in the cone used for this penetration test - auto cone_tip = DEMSim.AddWavefrontMeshObject(GetDEMEDataFile("mesh/cone.obj"), mat_type_cone); + // Load in the cone used for this penetration test. The penetrating tip is an analytical cone segment so the contact + // normal is exact and does not depend on the resolution of a triangle mesh. + float tip_height = std::sqrt(3.); + const float cone_radius = cone_diameter / 2; + const float cone_height = cone_radius * tip_height; + auto cone_tip = DEMSim.AddExternalObject(); + cone_tip->AddConeSegment(make_float3(0, 0, -3.0 / 4.0 * cone_height), make_float3(0, 0, 1), + 1.0 / tip_height, 0.0, cone_height, mat_type_cone, ENTITY_NORMAL_OUTWARD); auto cone_body = DEMSim.AddWavefrontMeshObject(GetDEMEDataFile("mesh/cyl_r1_h2.obj"), mat_type_cone); - std::cout << "Total num of triangles: " << cone_tip->GetNumTriangles() + cone_body->GetNumTriangles() << std::endl; + std::cout << "Total num of mesh triangles: " << cone_body->GetNumTriangles() << std::endl; - // The initial cone mesh has base radius 1, and height 1. Let's stretch it a bit so it has a 60deg tip, instead of - // 90deg. - float tip_height = std::sqrt(3.); - cone_tip->Scale(make_float3(1, 1, tip_height)); - // Then set mass properties - float cone_mass = 7.8e3 * tip_height / 3 * math_PI; + // Then set mass properties for the analytical cone segment. Its owner frame is placed at the solid cone centroid, + // so the local apex is at -3H/4 and the base rim at H/4. + float cone_mass = 7.8e3 * math_PI * cone_radius * cone_radius * cone_height / 3; cone_tip->SetMass(cone_mass); - // You can checkout https://en.wikipedia.org/wiki/List_of_moments_of_inertia - cone_tip->SetMOI(make_float3(cone_mass * (3. / 20. + 3. / 80. * tip_height * tip_height), - cone_mass * (3. / 20. + 3. / 80. * tip_height * tip_height), 3 * cone_mass / 10)); - // This cone mesh has its tip at the origin. And, float4 quaternion pattern is (x, y, z, w). - cone_tip->InformCentroidPrincipal(make_float3(0, 0, 3. / 4. * tip_height), make_float4(0, 0, 0, 1)); - // Note the scale method will scale mass and MOI automatically. But this only goes for the case you scale xyz all - // together; otherwise, the MOI scaling will not be accurate and you should manually reset them. - cone_tip->Scale(cone_diameter / 2); + cone_tip->SetMOI(make_float3(cone_mass * (3. / 20. * cone_radius * cone_radius + + 3. / 80. * cone_height * cone_height), + cone_mass * (3. / 20. * cone_radius * cone_radius + + 3. / 80. * cone_height * cone_height), + 3 * cone_mass * cone_radius * cone_radius / 10)); cone_tip->SetFamily(2); // The define the body that is connected to the tip @@ -214,11 +214,11 @@ int main() { // Put the cone in place double starting_height = terrain_max_z + 0.03; // Its initial position should be right above the cone tip... - body_tracker->SetPos(make_float3(0, 0, 0.5 + (cone_diameter / 2 / 4 * tip_height) + starting_height)); + body_tracker->SetPos(make_float3(0, 0, 0.5 + cone_height / 4 + starting_height)); // Note that position of objects is always the location of their centroid tip_tracker->SetPos(make_float3(0, 0, starting_height)); // The tip location, used to measure penetration length - double tip_z = -cone_diameter / 2 * 3 / 4 * tip_height + starting_height; + double tip_z = -3.0 / 4.0 * cone_height + starting_height; // Enable cone DEMSim.ChangeFamily(2, 1); @@ -235,8 +235,7 @@ int main() { for (float t = 0; t < sim_end; t += frame_time) { // float terrain_max_z = max_z_finder->GetValue(); float3 forces = tip_tracker->ContactAcc(); - // Note cone_mass is not the true mass, b/c we scaled the the cone tip! So we use true mass by using - // cone_tip->mass. + // ContactAcc returns acceleration for the analytical cone owner, so multiply by its assigned mass. forces *= cone_tip->mass; float pressure = std::abs(forces.z) / cone_surf_area; if (pressure > 1e-4 && !hit_terrain) { @@ -271,4 +270,4 @@ int main() { std::cout << "ConePenetration demo exiting..." << std::endl; return 0; -} \ No newline at end of file +} diff --git a/src/kernel/DEMHelperKernels.cuh b/src/kernel/DEMHelperKernels.cuh index d6a2729a..a73beef9 100644 --- a/src/kernel/DEMHelperKernels.cuh +++ b/src/kernel/DEMHelperKernels.cuh @@ -454,6 +454,16 @@ inline __host__ __device__ void matProxy2ContactParam(T1& E_eff, E_eff = (T1)1 / invE; } +inline __host__ __device__ float3 anyPerpendicularUnitVector(const float3& axis) { + float3 reference = (fabsf(axis.x) < 0.5f) ? make_float3(1, 0, 0) : make_float3(0, 1, 0); + float3 perpendicular = cross(axis, reference); + if (length(perpendicular) < DEME_TINY_FLOAT) { + reference = make_float3(0, 0, 1); + perpendicular = cross(axis, reference); + } + return normalize(perpendicular); +} + // Check whether a sphere and an analytical boundary are in contact, and gives overlap depth, contact point and contact // normal. template @@ -515,6 +525,64 @@ inline __host__ __device__ deme::contact_t checkSphereEntityOverlap(const T1& A, } return contactType; } + case (deme::ANAL_OBJ_TYPE_CONE_INF): + case (deme::ANAL_OBJ_TYPE_CONE): { + const float cone_slope = size1B; + const T3 min_h = (T3)size2B; + const T3 max_h = (T3)size3B; + const T1 tip2sph = A - B; + const T1 cone_axis = to_real3(dirB); + const T3 axial_dist = dot(tip2sph, cone_axis); + const T1 radial_vec = tip2sph - cone_axis * axial_dist; + const T3 radial_dist = length(radial_vec); + const float3 radial_dir = (radial_dist >= (T3)DEME_TINY_FLOAT) + ? to_real3(radial_vec / radial_dist) + : anyPerpendicularUnitVector(dirB); + + // The exact closest point on the conical side is found by projecting the sphere center onto the 2-D + // generator line (h, rho) = (h, slope * h). Clamping that axial coordinate gives the apex/rim feature for + // bounded cone segments without introducing mesh facets. + const T3 inv_slope_metric = (T3)1.0 / ((T3)1.0 + (T3)cone_slope * (T3)cone_slope); + T3 closest_h = (axial_dist + (T3)cone_slope * radial_dist) * inv_slope_metric; + bool closest_is_edge = false; + if (closest_h < min_h) { + closest_h = min_h; + closest_is_edge = true; + } else if (closest_h > max_h) { + closest_h = max_h; + closest_is_edge = true; + } + + if (closest_is_edge) { + const T1 closest = B + cone_axis * closest_h + + to_real3(radial_dir) * ((T3)cone_slope * closest_h); + const T1 feature2sph = A - closest; + const T3 dist_to_feature = length(feature2sph); + overlapDepth = (T3)radA + (T3)beta4Entity - dist_to_feature; + contactType = (overlapDepth < (T3)0.0) ? deme::NOT_A_CONTACT : deme::SPHERE_CONE_CONTACT; + + if (dist_to_feature >= (T3)DEME_TINY_FLOAT) { + cntNormal = to_real3(feature2sph / dist_to_feature); + } else { + const T3 side_normal_len = sqrt((T3)1.0 + (T3)cone_slope * (T3)cone_slope); + cntNormal = normal_sign * (cone_slope * dirB - radial_dir) / side_normal_len; + } + CP = A - to_real3(cntNormal * ((T3)radA - overlapDepth / (T3)2.0)); + return contactType; + } + + // Directional cone contact mirrors the cylinder convention. With inward normals, particles are expected + // where radial_dist <= slope * axial_dist; with outward normals, they are expected outside that side. + const T3 side_normal_len = sqrt((T3)1.0 + (T3)cone_slope * (T3)cone_slope); + const T3 signed_gap = + (T3)normal_sign * ((T3)cone_slope * axial_dist - radial_dist) / side_normal_len; + overlapDepth = (T3)radA + (T3)beta4Entity - signed_gap; + contactType = (overlapDepth < (T3)0.0) ? deme::NOT_A_CONTACT : deme::SPHERE_CONE_CONTACT; + + cntNormal = normal_sign * (cone_slope * dirB - radial_dir) / side_normal_len; + CP = A - to_real3(cntNormal * ((T3)radA - overlapDepth / (T3)2.0)); + return contactType; + } default: return deme::NOT_A_CONTACT; }