Hi Ruochun:

The script is attached below.

Thank you,
Weigang

On Thursday, June 20, 2024 at 7:27:35 PM UTC+8 Ruochun Zhang wrote:

> I assume you are doing it via a custom force model. In that case, this 
> problem is very similar to the exact demo we provided (which runs). You 
> have to tell us what you changed in order for us to give suggestions.
>
> Thank you,
> Ruochun
>
> On Thursday, June 20, 2024 at 6:55:10 PM UTC+8 [email protected] wrote:
>
>> Hi Ruochun,
>>
>> I am currently modelling the impact of a rock block against a plate. 
>> However, as shown below, the spheres penetrate through the plate without 
>> any contact with the plate. Could you give some suggestions?
>> [image: 微信图片_20240620185350.jpg]
>>
>> Thank you,
>> Weigang
>>
>> On Thursday, June 20, 2024 at 5:02:38 PM UTC+8 Ruochun Zhang wrote:
>>
>>> Hi Weigang,
>>>
>>> Yes. And that's exactly what we did in the fracture demo to begin with, 
>>> is it not?
>>>
>>> Thank you,
>>> Ruochun
>>>
>>> On Thursday, June 20, 2024 at 3:07:07 PM UTC+8 [email protected] 
>>> wrote:
>>>
>>>> Hi Ruochun,
>>>>
>>>> Do you mean that we can add a if statements in the 
>>>> ForceModelWithFracture.cu to  differentiate types of contacts?
>>>>
>>>> Thank you,
>>>> Weigang
>>>>
>>>>
>>>>
>>>> On Thursday, June 13, 2024 at 8:30:45 PM UTC+8 Ruochun Zhang wrote:
>>>>
>>>>> Currently, DEM-Engine custom model functionality permits only one 
>>>>> script, meaning that in principal, you are writing a unified force model 
>>>>> that goes for sphere-sphere, sphere-boundary, *and *sphere-mesh. Just 
>>>>> think of the sphere-compressor contact as a special case, which is the 
>>>>> contact between a normal sphere and an infinitely-large sphere.
>>>>>
>>>>> About enforcing different models for different objects, you can 
>>>>> associate different objects with respective family numbers, then use if 
>>>>> statements to differentiate types of contacts. I talked about this in one 
>>>>> of the earlier replies in this thread, please go back and have a look.
>>>>>
>>>>> Thank you,
>>>>> Ruochun
>>>>>
>>>>>
>>>>> On Thursday, June 13, 2024 at 1:31:37 PM UTC+8 [email protected] 
>>>>> wrote:
>>>>>
>>>>>> Hi Ruochun,
>>>>>>
>>>>>> Thank you for your reply. In fact, I am puzzled that there is not any 
>>>>>> sentence in the  Fracture-Box demo to declare the contact model between 
>>>>>> spheres and the compressing plate. So, how to declare the contact model, 
>>>>>> especially for the case in which there has been another contact model, 
>>>>>> such 
>>>>>> as the  ForceModelWithFracture.cu?
>>>>>>
>>>>>> Thank you,
>>>>>> Weigang
>>>>>>
>>>>>> On Thursday, June 13, 2024 at 12:03:52 AM UTC+8 Ruochun Zhang wrote:
>>>>>>
>>>>>>> Hi Weigang,
>>>>>>>
>>>>>>> Note that in general, if you use a custom force model, then that 
>>>>>>> model alone is used, there is no fallback so I wouldn't call that 
>>>>>>> anything 
>>>>>>> has "defaulted" to Hertz–Mindlin. Specific to the Fracture-Box demo, 
>>>>>>> it's 
>>>>>>> just that the force model (ForceModelWithFracture.cu) is written in a 
>>>>>>> way 
>>>>>>> that when the bond between two particles breaks, the "else" branch of 
>>>>>>> that 
>>>>>>> big if statement is executed, and something similar to Hertz–Mindlin is 
>>>>>>> enforced.
>>>>>>>
>>>>>>> I said it's similar since it's not a pure spring–damper model now, 
>>>>>>> it's a spring–damper model with the resting location (zero-force 
>>>>>>> location) 
>>>>>>> being the penetration depth the moment when the bond broke. This 
>>>>>>> implementation avoids a potential problem with a pure spring–damper 
>>>>>>> model: 
>>>>>>> Without the adjusted resting location, a newly broken-away pair of 
>>>>>>> particles would see the bond suddenly gone but there were still 
>>>>>>> left-over 
>>>>>>> penetration, therefore immediately pushing each other away with high 
>>>>>>> velocity, resulting in non-physical results. Shout-out to Bona for 
>>>>>>> implementing this.
>>>>>>>
>>>>>>> Thank you!
>>>>>>> Ruochun
>>>>>>>
>>>>>>>
>>>>>>> On Wednesday, June 12, 2024 at 8:33:55 PM UTC+8 [email protected] 
>>>>>>> wrote:
>>>>>>>
>>>>>>>> Hi Ruochun,
>>>>>>>>
>>>>>>>> Thank you for your great work. The demo now is runs well. I have a 
>>>>>>>> question about the Fracture-Box demo. After the bond between two 
>>>>>>>> particles 
>>>>>>>> breaks, is it right that the contact model between them change 
>>>>>>>> defaultly 
>>>>>>>> from ForceModelWithFracture to Hertz–Mindlin model?
>>>>>>>>
>>>>>>>> Thank you,
>>>>>>>> Weigang
>>>>>>>>
>>>>>>>> On Friday, May 31, 2024 at 8:05:24 PM UTC+8 Ruochun Zhang wrote:
>>>>>>>>
>>>>>>>>> Hi Weigang, 
>>>>>>>>>
>>>>>>>>> The problem appears to be fixed now. Please pull the newest main 
>>>>>>>>> branch and try the fracture demo again. I suggest you do a clean 
>>>>>>>>> rebuild 
>>>>>>>>> from an empty directory, because only a force model file is changed 
>>>>>>>>> and if 
>>>>>>>>> you just "make", cmake probably thinks nothing needs to be done. Or 
>>>>>>>>> you can 
>>>>>>>>> simply manually copy the modified cu file over to the appropriate 
>>>>>>>>> location 
>>>>>>>>> in the build folder.
>>>>>>>>>
>>>>>>>>> The problem was that we left an uninitialized variable in the 
>>>>>>>>> force model file. On data center cards, the compiler zeros it out 
>>>>>>>>> automatically, but on consumer cards it does not. We therefore missed 
>>>>>>>>> it 
>>>>>>>>> initially. Should be all good now.
>>>>>>>>>
>>>>>>>>> Thank you, 
>>>>>>>>> Ruochun
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Tue, May 28, 2024, 11:00 PM Ruochun Zhang <[email protected]> 
>>>>>>>>> wrote:
>>>>>>>>>
>>>>>>>>>> Hi Weigang,
>>>>>>>>>>
>>>>>>>>>> It's not like recommending a card; rather, if you happen to have 
>>>>>>>>>> a data center card available (A100, A5000 etc.), you can circumvent 
>>>>>>>>>> this 
>>>>>>>>>> problem for now by using them. Otherwise, let's hope we find a 
>>>>>>>>>> solution 
>>>>>>>>>> very soon, or at least find a workaround. We'll keep you posted.
>>>>>>>>>>
>>>>>>>>>> Thank you,
>>>>>>>>>> Ruochun
>>>>>>>>>>
>>>>>>>>>> On Monday, May 27, 2024 at 7:27:10 PM UTC+8 [email protected] 
>>>>>>>>>> wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi Ruochun,
>>>>>>>>>>>
>>>>>>>>>>> My card is indeed a gaming card NVIDIA GeForce RTX 2070 Super. 
>>>>>>>>>>> Could you recommend a card?
>>>>>>>>>>>
>>>>>>>>>>> Thanks,
>>>>>>>>>>> Weigang
>>>>>>>>>>>
>>>>>>>>>>> On Saturday, May 25, 2024 at 10:44:17 PM UTC+8 Ruochun Zhang 
>>>>>>>>>>> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> Hi Weigang,
>>>>>>>>>>>>
>>>>>>>>>>>> I would like to answer Question 2 first. You can use variables 
>>>>>>>>>>>> *AOwnerFamily 
>>>>>>>>>>>> *and *BOwnerFamily *to refer to clump (or owner) A's and clump 
>>>>>>>>>>>> B's family numbers directly in the force model. So you can write 
>>>>>>>>>>>> *if 
>>>>>>>>>>>> *statements that treat each case accordingly. This is 
>>>>>>>>>>>> mentioned in the DEME paper: You can check out Table 2.
>>>>>>>>>>>>
>>>>>>>>>>>> For Question 1, the situation is more complicated. We are 
>>>>>>>>>>>> actually able to reproduce the problem you encountered (simulation 
>>>>>>>>>>>> freezing), but on consumer/gaming cards only. It runs perfectly on 
>>>>>>>>>>>> data 
>>>>>>>>>>>> center cards and that's what we used for generating this 
>>>>>>>>>>>> experiment, 
>>>>>>>>>>>> therefore it did not appear a problem for us. It's likely due to 
>>>>>>>>>>>> different 
>>>>>>>>>>>> GPU architectures running the same code differently, we'll try to 
>>>>>>>>>>>> find a 
>>>>>>>>>>>> solution and let you know in the following days, so stay tuned. By 
>>>>>>>>>>>> the way, 
>>>>>>>>>>>> you are using a gaming card to run the simulations, correct?
>>>>>>>>>>>>
>>>>>>>>>>>> Thank you,
>>>>>>>>>>>> Ruochun
>>>>>>>>>>>>
>>>>>>>>>>>> On Friday, May 24, 2024 at 8:47:00 AM UTC+8 [email protected] 
>>>>>>>>>>>> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>> Thank you for your reply! I have two questions about the demo 
>>>>>>>>>>>>> "DEMdemo-Fracture-Box".
>>>>>>>>>>>>>
>>>>>>>>>>>>> (1) After i have run this demo for one day, there is still 
>>>>>>>>>>>>> nothing outputed. Is it fact that the DEME will run a long time 
>>>>>>>>>>>>> once a 
>>>>>>>>>>>>> fracture model is used?
>>>>>>>>>>>>>
>>>>>>>>>>>>> (2) How to tell DEME the contact force models between 
>>>>>>>>>>>>> different familys? I did not see any snippet which undertakes 
>>>>>>>>>>>>> this work, 
>>>>>>>>>>>>> although the demo has at least two kinds of contact force model. 
>>>>>>>>>>>>> The 
>>>>>>>>>>>>> question i really want to ask is that, if we have more than three 
>>>>>>>>>>>>> familys 
>>>>>>>>>>>>> in our model and there are more than two familys consist of 
>>>>>>>>>>>>> spherical 
>>>>>>>>>>>>> particles, how to tell DEME the  contact force model between 
>>>>>>>>>>>>> different 
>>>>>>>>>>>>> familys, and tell DEME the contact force model between the 
>>>>>>>>>>>>> particles of a 
>>>>>>>>>>>>> family in the same time.
>>>>>>>>>>>>>
>>>>>>>>>>>>> Thank you,
>>>>>>>>>>>>> Weigang
>>>>>>>>>>>>>
>>>>>>>>>>>>> On Thursday, May 23, 2024 at 10:10:14 PM UTC+8 Ruochun Zhang 
>>>>>>>>>>>>> wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> No. That's because the *sphere components* in a clump are 
>>>>>>>>>>>>>> just shape representations, and they have no physics, only 
>>>>>>>>>>>>>> material 
>>>>>>>>>>>>>> properties. If they did have interaction with each other then 
>>>>>>>>>>>>>> every clump 
>>>>>>>>>>>>>> would've immediately exploded upon simulation starting.
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Thank you,
>>>>>>>>>>>>>> Ruochun
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Thursday, May 23, 2024 at 12:16:32 PM UTC+8 
>>>>>>>>>>>>>> [email protected] wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Hi Ruochun,
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Thank you for your reply! I have a question about the clump. 
>>>>>>>>>>>>>>> Do the particles in a clump interact with each other via any 
>>>>>>>>>>>>>>> force model?
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Thank you,
>>>>>>>>>>>>>>> Weigang
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> On Monday, May 20, 2024 at 10:26:03 PM UTC+8 Ruochun Zhang 
>>>>>>>>>>>>>>> wrote:
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Hi Weigang,
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Yes. Again,  DEME only cares about the position and radius 
>>>>>>>>>>>>>>>> of the initial spheres, so you can supply it however you want. 
>>>>>>>>>>>>>>>> If you have 
>>>>>>>>>>>>>>>> a file that records this information using your own format, 
>>>>>>>>>>>>>>>> then the 
>>>>>>>>>>>>>>>> easiest thing is to have your own C++ or Python script that 
>>>>>>>>>>>>>>>> reads it into 
>>>>>>>>>>>>>>>> arrays and then feeds them to the solver. If your file 
>>>>>>>>>>>>>>>> shares the format with DEME's standard clump output file, or 
>>>>>>>>>>>>>>>> your file is 
>>>>>>>>>>>>>>>> simply the output of *WriteClumpFile*, then you can 
>>>>>>>>>>>>>>>> conveniently load clump positions and orientations by 
>>>>>>>>>>>>>>>> *ReadClumpXyzFromCsv 
>>>>>>>>>>>>>>>> *and *ReadClumpQuatFromCsv *(examples are in GRCPrep_Part1 
>>>>>>>>>>>>>>>> and 2). Note that currently, sphere radius cannot be read from 
>>>>>>>>>>>>>>>> standard 
>>>>>>>>>>>>>>>> clump output files since it's part of the clump template 
>>>>>>>>>>>>>>>> information thus 
>>>>>>>>>>>>>>>> not in the clump output.
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Thank you,
>>>>>>>>>>>>>>>> Ruochun
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> On Monday, May 20, 2024 at 12:56:56 PM UTC+8 
>>>>>>>>>>>>>>>> [email protected] wrote:
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Hi Everyone,
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> I would like to generate a block which consists of a 
>>>>>>>>>>>>>>>>> assembly spheres. Is it possible to generate the block by 
>>>>>>>>>>>>>>>>> other code, and 
>>>>>>>>>>>>>>>>> then load it into DEM-Engine via a file containg the 
>>>>>>>>>>>>>>>>> informations (such as 
>>>>>>>>>>>>>>>>> position and radius) of the spheres?
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> Thank you,
>>>>>>>>>>>>>>>>> Weigang
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> -- 
>>>>>>>>>> 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/1910af21-d1e3-43dd-a779-0c229e6150a6n%40googlegroups.com
>>>>>>>>>>  
>>>>>>>>>> <https://groups.google.com/d/msgid/projectchrono/1910af21-d1e3-43dd-a779-0c229e6150a6n%40googlegroups.com?utm_medium=email&utm_source=footer>
>>>>>>>>>> .
>>>>>>>>>>
>>>>>>>>>

-- 
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/f0d53f93-3d0d-439e-b1c8-4c02030144a5n%40googlegroups.com.
//  Copyright (c) 2021, SBEL GPU Development Team
//  Copyright (c) 2021, University of Wisconsin - Madison
//
//	SPDX-License-Identifier: BSD-3-Clause

// =============================================================================
// Fracture
// =============================================================================

#include <DEM/API.h>
#include <DEM/HostSideHelpers.hpp>
#include <DEM/utils/Samplers.hpp>

#include <chrono>
#include <cmath>
#include <cstdio>
#include <filesystem>
#include <iostream>
#include <fstream>

using namespace deme;

const double math_PI = 3.1415927;

int main() {
    DEMSolver DEMSim;
    DEMSim.SetVerbosity(INFO);
    DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);
    DEMSim.SetOutputContent(OUTPUT_CONTENT::VEL);
    DEMSim.SetMeshOutputFormat(MESH_FORMAT::VTK);
    DEMSim.SetContactOutputContent(OWNER | FORCE | CNT_WILDCARD);

    // This demo could lead to large numbers of per-sphere contacts, so to be safe...
    DEMSim.SetErrorOutAvgContacts(500);

    //  E, nu, CoR, mu, Crr...
    auto mat_type_container =
        DEMSim.LoadMaterial({{"E", 1e8}, {"nu", 0.3}, {"CoR", 0.7}, {"mu", 0.80}, {"Crr", 0.10}});
    auto mat_type_particle = DEMSim.LoadMaterial({{"E", 1e8}, {"nu", 0.20}, {"CoR", 0.5}, {"mu", 0.5}, {"Crr", 0.05}});
    // If you don't have this line, then values will take average between 2 materials, when they are in contact
    DEMSim.SetMaterialPropertyPair("CoR", mat_type_container, mat_type_particle, 0.2);
    DEMSim.SetMaterialPropertyPair("Crr", mat_type_container, mat_type_particle, 0.5);
    DEMSim.SetMaterialPropertyPair("mu", mat_type_container, mat_type_particle, 0.5);
    // We can specify the force model using a file.
    auto my_force_model = DEMSim.ReadContactForceModel("ForceModelWithFractureModel.cu");

    // Those following lines are needed. We must let the solver know that those var names are history variable etc.
    my_force_model->SetMustHaveMatProp({"E", "nu", "CoR", "mu", "Crr"});
    my_force_model->SetMustPairwiseMatProp({"CoR", "mu", "Crr"});
    // Pay attention to the extra per-contact wildcard `unbroken' here.
    my_force_model->SetPerContactWildcards(
        {"delta_time", "delta_tan_x", "delta_tan_y", "delta_tan_z", "unbroken", "initialLength"});

    float world_size = 1.5;
    float container_diameter = 0.05;
    float terrain_density = 2.80e3;
    float sphere_rad = 0.0020;

    float step_size = 2e-7;
    float fact_radius = 0.9;

    DEMSim.InstructBoxDomainDimension(world_size, world_size, world_size);
    // No need to add simulation `world' boundaries, b/c we'll add a cylinderical container manually
    DEMSim.InstructBoxDomainBoundingBC("all", mat_type_container);
    // DEMSim.SetInitBinSize(sphere_rad * 5);
    //  Now add a cylinderical boundary along with a bottom plane
    double bottom = -0;
    double top = 0.10;

    // auto walls = DEMSim.AddExternalObject();
    // std::string shake_pattern_xz = " 0.00 * sin( 500 * 2 * deme::PI * t)";
    // walls->AddPlane(make_float3(0, 0, bottom), make_float3(0, 0, 1), mat_type_container);
    // walls->AddPlane(make_float3(0, 0, world_size / 2. - world_size / 20.), make_float3(0, 0, -1),
    // mat_type_container); walls->SetFamily(2);

    auto walls = DEMSim.AddWavefrontMeshObject("../data/mesh/funnel_left.obj", mat_type_container);
    float3 move = make_float3(0.05, 0.00, 0 - 2 * sphere_rad);  // z
    float4 rot = make_float4(0.7071, 0, 0, 0.7071);
    walls->Scale(make_float3(0.8, 0.1, 2.0));
    walls->Move(move, rot);
    walls->SetFamily(2);
    DEMSim.SetFamilyFixed(2);

    auto cylinder = DEMSim.AddExternalObject();
    cylinder->AddCylinder(make_float3(0), make_float3(0, 0, 1), 1.8 * container_diameter / 2., mat_type_container, 0);
    cylinder->SetFamily(10);
    DEMSim.SetFamilyFixed(10);
    // DEMSim.SetFamilyPrescribedLinVel(10, shake_pattern_xz, shake_pattern_xz, shake_pattern_xz);

    //auto plate = DEMSim.AddWavefrontMeshObject("../data/mesh/funnel_left.obj", mat_type_container);
    //move = make_float3(0.05, 0.00, top + 2.0 * sphere_rad);  // z
    //rot = make_float4(0.7071, 0, 0, 0.7071);
    //plate->Scale(make_float3(0.8, 0.1, 2.0));
    //plate->Move(move, rot);
    //plate->SetFamily(20);
    //DEMSim.SetFamilyFixed(20);

    // Track it

    auto plate_tracker = DEMSim.Track(walls);

    //DEMSim.SetFamilyPrescribedLinVel(22, "0", "0", to_string_with_precision(-0.0050));

    // auto total_mass_finder = DEMSim.CreateInspector("clump_mass");
    // auto max_v_finder = DEMSim.CreateInspector("clump_max_absv");

    // Define the terrain particle templates

    // Calculate its mass and MOI

    float sphere_vol = 4. / 3. * math_PI * sphere_rad * sphere_rad * sphere_rad;
    float mass = terrain_density * sphere_vol;
    // Then load it to system
    std::shared_ptr<DEMClumpTemplate> my_template = DEMSim.LoadSphereType(mass, sphere_rad, mat_type_particle);

    // Sampler to sample
    HCPSampler sampler(2.0 * sphere_rad);

    float fill_height = top;
    float3 fill_center = make_float3(0, 0, bottom + fill_height / 2);
    const float fill_radius = container_diameter / 2. - sphere_rad * 0.;
    float3 fill = make_float3(fill_radius, fill_radius, bottom + fill_height / 2);
    auto input_xyz = sampler.SampleBox(fill_center, fill);
    // auto input_xyz = sampler(fill_center, fill_radius, fill_height / 2 - sphere_rad * 0.);
    auto particles = DEMSim.AddClumps(my_template, input_xyz);
    particles->SetFamily(1);
    std::cout << "Total num of particles: " << particles->GetNumClumps() << std::endl;
    
    DEMSim.SetFamilyPrescribedLinVel(1, "0", "0", to_string_with_precision(-1.0));

    std::filesystem::path out_dir = std::filesystem::current_path();
    std::string nameOutFolder = "R" + std::to_string(sphere_rad) + "_Int" + std::to_string(fact_radius) + "";
    out_dir += "/BlockImpact_" + nameOutFolder;
    remove_all(out_dir);
    create_directory(out_dir);

    // Some inspectors

    auto max_z_finder = DEMSim.CreateInspector("clump_max_z");
    auto min_z_finder = DEMSim.CreateInspector("clump_min_z");

    DEMSim.SetFamilyExtraMargin(1, fact_radius * sphere_rad);

    DEMSim.SetInitTimeStep(step_size);
    DEMSim.SetGravitationalAcceleration(make_float3(0, 0.00, 1 * -9.81));
    DEMSim.Initialize();
    // DEMSim.DisableContactBetweenFamilies(20, 1);
    std::cout << "Initial number of contacts: " << DEMSim.GetNumContacts() << std::endl;

    float sim_end = 5;
    unsigned int fps = 1000;
    unsigned int datafps = 1000;
    unsigned int modfpsGeo = datafps / fps;
    float frame_time = 1.0 / datafps;
    std::cout << "Output at " << fps << " FPS" << std::endl;
    unsigned int out_steps = (unsigned int)(1.0 / (datafps * step_size));
    unsigned int frame_count = 0;
    unsigned int step_count = 0;

    bool status_1 = true;
    bool status_2 = true;

    // DEMSim.DisableContactBetweenFamilies(10, 1);

    double L0;
    double stress;
    std::string nameOutFile = "data_R" + std::to_string(sphere_rad) + "_Int" + std::to_string(fact_radius) + ".csv";
    std::ofstream csvFile(nameOutFile);

    DEMSim.SetFamilyContactWildcardValueAll(1, "initialLength", 0.0);
    // DEMSim.SetFamilyContactWildcardValueAll(1, "damage", 0.0);
    DEMSim.SetFamilyContactWildcardValueAll(1, "unbroken", 0.0);

    // Simulation loop
    for (float t = 0; t < sim_end; t += frame_time) {
        // DEMSim.ShowThreadCollaborationStats();

        std::cout << "Initial number of contacts: " << DEMSim.GetNumContacts() << std::endl;


        if (frame_count % 1 == 0) {
            float3 forces = plate_tracker->ContactAcc();
            double L = max_z_finder->GetValue() - min_z_finder->GetValue() + 2 * sphere_rad;
            std::cout << "Time: " << t << std::endl;
            std::cout << "Force [N]: " << stress << std::endl;

            if (frame_count % modfpsGeo == 0) {
                char filename[200];
                char meshname[200];
                char cnt_filename[200];

                std::cout << "Outputting frame: " << frame_count / modfpsGeo << std::endl;
                sprintf(filename, "%s/DEMdemo_output_%04d.csv", out_dir.c_str(), frame_count / modfpsGeo);
                sprintf(meshname, "%s/DEMdemo_mesh_%04d.vtk", out_dir.c_str(), frame_count / modfpsGeo);
                sprintf(cnt_filename, "%s/DEMdemo_contact_%04d.csv", out_dir.c_str(), frame_count / modfpsGeo);

                DEMSim.WriteSphereFile(std::string(filename));
                DEMSim.WriteMeshFile(std::string(meshname));
                DEMSim.WriteContactFile(std::string(cnt_filename));
            }
        }

        DEMSim.DoDynamics(frame_time);
        frame_count++;
    }
    csvFile.close();
    DEMSim.ShowTimingStats();
    std::cout << "Fracture demo exiting..." << std::endl;
    return 0;
}

Reply via email to