Hi Wenxuan, Regarding debugging on Linux, I recommend seeking professional assistance. However, I find that VSCode's GUI for stepping through code is sufficient for my needs.
For enforcing motion on an entity step-by-step, use *DoDynamics(step_size)* or *DoStepDynamics()* instead of *DoDynamics(frame_duration)*. This way, you can count the number of steps and use if statements to control outputs. The RotatingDrum demo provides a good example of this. The performance impact should be minimal, as DEME is designed for fine-grain operations like these to enable more flexible co-simulation. The two GPU streams synchronize only with *DoDynamicsThenSync()*, not *DoDynamics()*. Also, transferring the explicit position and velocity of one or two entities to GPUs at each step introduces only a small and manageable overhead, too. Regarding your script, it doesn't fully implement my specified formula. I didn't use *fabs*, and the linear velocity should have a scaling factor of 5, which seems to be missing. Even with my specified formula, I can't guarantee its correctness (without knowing details like what the coordinate system is, and what the "expected results" is). But more realistically, it's that I cannot do your research for you. You see, this is becoming more of a consulting task than troubleshooting. I suggest starting simple, using only one of the implemented motions, verifying it works as expected, and then gradually adding components of the second one to approach your desired outcome. Thank you, Ruochun On Saturday, June 29, 2024 at 2:16:35 AM UTC+8 [email protected] wrote: > Hi, Ruochun > > I have tried it, but the geometric motion doesn't seem to achieve > the expected result. Can you help me troubleshoot the error? Ths mesh obj > file are attached. > > #include <core/ApiVersion.h> > #include <core/utils/ThreadManager.h> > #include <DEM/API.h> > #include <DEM/HostSideHelpers.hpp> > #include <DEM/utils/Samplers.hpp> > > #include <cstdio> > #include <chrono> > #include <filesystem> > > using namespace deme; > using namespace std::filesystem; > > int main() > { > DEMSolver DEMSim; > DEMSim.SetVerbosity(STEP_METRIC); > // For general use cases, you want to set the verbosity to INFO: It's also > a bit faster than STEP_METRIC. > // DEMSim.SetVerbosity(INFO); > DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV); > DEMSim.SetOutputContent(OUTPUT_CONTENT::ABSV); > DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK); > > // If you don't need individual force information, then this option makes > the solver run a bit faster. > DEMSim.SetNoForceRecord(); > > // E, nu, CoR, mu, Crr... > auto mat_type_mixer = DEMSim.LoadMaterial({{"E", 1e8}, {"nu", 0.3}, {"CoR", > 0.6}, {"mu", 0.5}, {"Crr", 0.0}}); > auto mat_type_granular = DEMSim.LoadMaterial({{"E", 1e8}, {"nu", 0.3}, { > "CoR", 0.6}, {"mu", 0.2}, {"Crr", 0.0}}); > // If you don't have this line, then mu between mixer material and > granular material will be 0.35 (average of the > // two). > DEMSim.SetMaterialPropertyPair("mu", mat_type_mixer, mat_type_granular, > 0.5); > > float step_size = 5e-5; > const double world_size = 10; > const float chamber_height = world_size / 3.; > const float fill_height = chamber_height; > const float chamber_bottom = -world_size / 2.; > const float fill_bottom = chamber_bottom + chamber_height; > > DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size); > DEMSim.InstructBoxDomainBoundingBC("all", mat_type_granular); > > > auto mixer = DEMSim.AddWavefrontMeshObject((GET_DATA_PATH() / > "mesh/Chute_Rotating.obj").string(), mat_type_mixer); > std::cout << "Total num of triangles: " << mixer->GetNumTriangles() << std > ::endl; > mixer->SetFamily(10); > mixer->InformCentroidPrincipal(make_float3(0, 0, 6.793), make_float4(0, 0, > 0, 1)); > > > double w1 = 3.14159; > double w2 = 1.0; > std::string angVel_x = "fabs(1.0 * cos(deme::PI * t))"; > std::string angVel_y = "fabs(1.0 * sin(deme::PI * t))"; > std::string angVel_z = "deme::PI"; > > > std::string vel_x = "(1.0 * cos(1.0 * t) * sin(deme::PI * t) + sin(1.0 * > t) * deme::PI * cos(deme::PI * t))"; > std::string vel_y = "(-1.0 * cos(1.0 * t) * cos(deme::PI * t) + sin(1.0 * > t) * deme::PI * sin(deme::PI * t))"; > std::string vel_z = "0"; > > DEMSim.SetFamilyPrescribedAngVel(10, angVel_x, angVel_y, angVel_z); > DEMSim.SetFamilyPrescribedLinVel(10, vel_x, vel_y, vel_z); > > > > > float granular_rad = 0.005; > // Calculate its mass and MOI > float mass = 2.6e3 * 5.5886717; // in kg or g > float3 MOI = make_float3(2.928, 2.6029, 3.9908) * 2.6e3; > std::shared_ptr<DEMClumpTemplate> template_granular = DEMSim.LoadClumpType > (mass, MOI, GetDEMEDataFile("clumps/3_clump.csv"), mat_type_granular); > template_granular->Scale(granular_rad); > > // Track the mixer > auto mixer_tracker = DEMSim.Track(mixer); > > > > DEMSim.SetInitTimeStep(step_size); > DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81)); > > // Initialize the simulation system > DEMSim.Initialize(); > > path out_dir = current_path(); > out_dir += "/DEMdemo_RotatingChute"; > create_directory(out_dir); > > float sim_end = 2.0; > unsigned int fps = 20; > float frame_time = 1.0 / fps; > > std::cout << "Output at " << fps << " FPS" << std::endl; > unsigned int currframe = 0; > > > // mixer_tracker->SetOriQ(make_float4(0, 0, 0, 1)); > for (float t = 0; t < sim_end; t += frame_time) // 0.05迭代一次 > { > > std::cout << "Position: " << mixer_tracker->GetPos()[0] << " " << > mixer_tracker->GetPos()[1] << " " << mixer_tracker->GetPos()[2] << std:: > endl; > std::vector<float> Qua = mixer_tracker->GetOriQ(); // x, y, z, w > // std::cout << Qua[0] << " " << Qua[1] << " " << Qua[2] << " " << Qua[3] > << std::endl; > > > std::cout << "Frame: " << currframe << std::endl; > char filename[200], meshfilename[200], cnt_filename[200]; > sprintf(filename, "%s/DEMdemo_output_%04d.csv", out_dir.c_str(), currframe > ); > sprintf(meshfilename, "%s/DEMdemo_mesh_%04d.vtk", out_dir.c_str(), > currframe); > sprintf(cnt_filename, "%s/Contact_pairs_%04d.csv", out_dir.c_str(), > currframe++); > > DEMSim.WriteMeshFile(std::string(meshfilename)); > DEMSim.DoDynamics(frame_time); > > } > > std::cout << "DEMdemo_RotatingChute simulating complete!" << std::endl; > return 0; > } > > 在2024年6月29日星期六 UTC+8 00:33:30<simboys> 写道： > >> Hi, Ruochun >> >> Thank you for your reply, I wiil try it. By the way, can you share >> me the method how to debug the code step-by-step in linux? In addition, if >> I use the DEMTracker to set the location and rotation of geometry at each >> time step, will this affect computational efficiency? And how to set the >> location of geometry at each time step? From the >> demo, DoDynamics(frameTime) is many times longer than timestep, and if we >> set the frameTime by timestep, it will result in a huge data file. >> >> Best regards, >> Wenxuan >> >> 在2024年6月27日星期四 UTC+8 23:06:58<Ruochun Zhang> 写道： >> >>> Hi Wenxuan, >>> >>> 1. When you load a mesh or an analytical object, it's by default in >>> family 255, a reserved and automatically fixed family. You probably changed >>> your mesh into a family that does not have any motion prescription, then in >>> that case it follows the "normal physics" like all the particles you load >>> in, including falling under gravity. If you wish to fix that family, just >>> use *SetFamilyFixed*. >>> >>> 2. This is more complicated. First, we have to understand that in DEME, >>> strictly speaking, objects work with their principal frames of inertia >>> (centered at CoM, having a diagonal inertia matrix). This is a current >>> limitation due to how the inertia matrix-defining APIs are designed: you >>> can only input the principal moments of inertia. But this also means that >>> setting the frame correctly is not necessary if you completely prescribe >>> the motion of that object (a.k.a. this object does not respond to physics), >>> and you don't care about getting a torque reading. However, in your case, I >>> suggest you do so, since you'll prescribe a complex motion, and >>> standardizing the frame to use at least the CoM frame can save you some >>> head-scratching. I suggest that you make it such that in the mesh file you >>> load in, the (0,0,0) point is right at the CoM of this mesh; but if not, >>> you can adjust the CoM and principal axes with *InformCentroidPrincipal* >>> after loading in the mesh. >>> >>> After that is taken care of, you need to combine the two rotations into >>> a rotation about this object's own CoM and a linear translation of this >>> object's CoM (because family motion prescription allows for setting these >>> two motions only). Use your MBD knowledge to do that. I don't guarantee the >>> correctness of this, but some back-of-the-envelop calculation seems to put >>> the angular velocity to be (w2×cos(w1×t), w2×sin(w1×t), w1), and the linear >>> velocity to be 5×(w2×cos(w2×t)×sin(w1×t)+sin(w2×t)×w1×cos(w1×t), >>> -w2×cos(w2×t)×cos(w1×t)+sin(w2×t)×w1×sin(w1×t), 0). Again, please verify >>> yourself. After that, prescribe the object's family's motion using it. >>> >>> But if this is too complicated and you can write your code to calculate >>> the location and quaternion of this object at any given *t *(which for >>> some people is easier than mathematically deriving the analytical >>> expression for the linear and rotational motions), then you can set this >>> object fixed, then use the tracker methods *SetPos *and *SetOriQ *to >>> enforce the location and rotation of this object *at each time step*. >>> It should achieve the same assuming this object is not moving at a very >>> high velocity. >>> >>> Thank you, >>> Ruochun >>> >>> On Thursday, June 27, 2024 at 2:49:43 PM UTC+8 [email protected] wrote: >>> >>>> Hi, >>>> I want to achieve a composite rotational motion of a geometric >>>> body, as follows >>>> Motion 1: Rotate at a speed of w1 rad/s, with the rotation axis being >>>> the z-axis >>>> Motion 2: Rotate at a speed of w2 rad/s, with a rotation vector of [1 * >>>> cos (w1 * time), 1 * sin (w1 * time), 0], and an action point of [0, 0, >>>> 5], >>>> as shown in attached figure. >>>> [image: rotating_Chute.png] >>>> >>>> >>>> However, the following issues were encountered: >>>> 1. After using DEMTracker to track geometry, even if I don't add any >>>> motion, the geometry will continue to drop down. How can I fix the >>>> geometry? >>>> 2. How to achieve the composite rotational motion mentioned above? >>>> >>>> Best regards, >>>> Wenxuan Xu >>>> >>> -- You received this message because you are subscribed to the Google Groups "ProjectChrono" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/projectchrono/2d801bf8-978d-421a-97cb-d7982eda636fn%40googlegroups.com.
