Different rebound againt analytic-box and fixed triangle-mesh wall

[박선호]

I am testing DEME mesh-wall contact with pyDEME.

I created two single-sphere cases intended to differ only in the wall representation:

1. Analytic box:
   InstructBoxDomainDimension(...) +
   InstructBoxDomainBoundingBC("all", wall_material)

2. Fixed triangle mesh box:
   AddWavefrontMeshObject(...) using a closed, manifold, inward-facing OBJ
   with the same 50 x 50 x 50 mm internal dimensions.

Both cases use exactly the same particle, material-pair, solver, and initial
conditions:

• one sphere, radius = 1 mm
• density = 2600 kg/m^3
• initial position = [10, 20, 10] mm
• gravity = [0, -9.81, 0] m/s^2
• particle material:
  E = 1e8 Pa, nu = 0.30, CoR = 0.80, mu = 0.20, Crr = 0
• wall material:
  E = 1e8 Pa, nu = 0.30, CoR = 0.60, mu = 0.50, Crr = 0
• particle-wall pair:
  CoR = 0.50, mu = 0.35
• initial timestep = 1e-6 s
• CD update frequency = 40
• SetExpandSafetyAdder(2.0)
• SetCDNumStepsMaxDriftMultipleOfAvg(1.2)
• SetCDNumStepsMaxDriftAheadOfAvg(6)
• SetSortContactPairs(True)
• output interval = 0.02 s

The predicted first-impact speed is about 0.610557 m/s. Before impact, both
cases match exactly:

| Time | Analytic box | Fixed mesh box |
| 0.06 s | y = 2.34405 mm, speed = 0.588672 m/s | y = 2.34405 mm, speed = 0.588672 m/s |

However, after the first floor collision they differ substantially:

| Time | Analytic box | Fixed mesh box |
| 0.08 s | y = 5.00014 mm, speed = 0.140174 m/s | y = 10.70100 mm, speed = 0.460826 m/s |

The analytic case then settles normally. The fixed triangle-mesh case has a
much stronger rebound. With a larger 5e-6 s timestep, repeated mesh impacts
eventually increase velocity until the solver aborts at the 20 m/s velocity
limit. Reducing to 1e-6 s prevents the immediate instability but the rebound
difference remains.

[Ruochun Zhang]

This is an interesting question. The first thing I'd ask you to do is to share with us the "01_box_contact_comparison_1particle.py" script if possible. It appears that this is your main comparison script and you did not share it with us.

Then, without running the script, I think if you used the mesh without any scaling, aka a 50m × 50m × 50m box in the unit system you are using, then because its size is too huge compared to the 1mm sphere, the floating-point-stored mesh XYZ information incurs significant inaccuracy when the sphere--triangle penetration derivation is executed. That is to say, you won't get an accurate penetration value if it's a small sphere set on an incredibly huge triangle facet, from our mixed-precision algorithm design perspective, so you just have to avoid that.

If you scaled the mesh before using it, so it is about the same size as the 1mm particle, then it's more interesting. Though DEME2's facet--sphere-based penetration can induce some inaccuracy when the sphere is in contact with multiple facets, I don't think it would be this large. So I'd see the test script first.

Thank you,

Ruochun

[박선호]

Thank you. I did scale the OBJ from mm to m before using it:

    mesh_wall.Scale([1e-3, 1e-3, 1e-3])

So the mesh is a 0.05 x 0.05 x 0.05 m box, with a sphere radius of 0.001 m, not a 50 m box.

I attached the missing main script:
- 01_box_contact_comparison_1particle.py

The analytic and mesh cases share the same sphere, material pair, gravity, initial position, timestep, and CD settings. Only the wall definition changes.

The mesh is closed/manifold with inward normals: 0 boundary edges, 0 non-manifold edges, signed volume `-125000 mm^3` before scaling.

At `t=0.06 s`, both cases match: `v=0.588672 m/s`.
At `t=0.08 s`, after the first floor collision:
- analytic: `v=0.140174 m/s`
- mesh: `v=0.460826 m/s`

I also reduced the timestep to `0.1 us`; the mesh rebound decreased slightly but remained about `0.431 m/s`.

Could multiple coplanar floor facets be responsible, or is there a recommended mesh-contact setting for this closed-container case?

2026년 8월 13일 목요일 오전 12시 3분 0초 UTC+9에 Ruochun Zhang님이 작성:

[Ruochun Zhang]

I have found the cause. This is unfortunately a mixed-precision issue. The triangle facets in the simulation have ~0.001 edge length, and 32-bit floating-point-based normal/distance calculations are somewhat unstable in this region. The quickest way to see this effect is to change MM_TO_M to 1e-1 (enlarge the mesh 100 times, no change to the physics), then you should see the analytical and meshed results become the same. The more complete fix is to change your implicit unit system from SI to a mm-based or cm-based length scale; that'd also fix the problem. I attached the sphere-bouncing results in the meshed case, after the mm-based unit system change.

Frame 0: t=0.000s max_v=0.000000 avg_contacts=0.000000 update_freq=40.0
Frame 1: t=0.020s max_v=196.232773 avg_contacts=0.000000 update_freq=29.0
Frame 2: t=0.040s max_v=392.258362 avg_contacts=0.000000 update_freq=28.5
Frame 3: t=0.060s max_v=588.428345 avg_contacts=1.000000 update_freq=28.0
Frame 4: t=0.080s max_v=140.121323 avg_contacts=0.000000 update_freq=28.5
Frame 5: t=0.100s max_v=56.106407 avg_contacts=0.000000 update_freq=28.5
Frame 6: t=0.120s max_v=252.344238 avg_contacts=0.000000 update_freq=28.0
Frame 7: t=0.140s max_v=25.544058 avg_contacts=1.000000 update_freq=27.5
Frame 8: t=0.160s max_v=71.380043 avg_contacts=1.000000 update_freq=27.0
Frame 9: t=0.180s max_v=1.446599 avg_contacts=1.000000 update_freq=27.5
Frame 10: t=0.200s max_v=0.000390 avg_contacts=1.000000 update_freq=25.0
Frame 11: t=0.220s max_v=0.000390 avg_contacts=1.000000 update_freq=25.0
Frame 12: t=0.240s max_v=0.000390 avg_contacts=1.000000 update_freq=25.0

It seems the penetration derivation is more sensitive to tiny mesh edges than I originally thought. But DEME is a core solver and it's dimensionless, so you can always pick your own unit system. In summary, I suggest switching to cm or mm to circumvent this issue.

Also note that in the upcoming DEME version 3.0, the triangle contact precision selection is changed and as I was testing, it appears to offer more accuracy for ~0.001 edge length cases, though still not completely the same as the analytical boundary. I may follow up later in this thread to invite you to test the same problem using DEME3.

Thank you,

Ruochun

[박선호]

Hi Ruochun,

I tested the same 1 mm sphere in a 50 x 50 x 50 mm box using consistent mm- and cm-based unit systems.

The SI mesh case showed excessive rebound and eventual instability. After changing to mm/cm units, the instability disappeared, but the mesh rebound still differs substantially from the analytical wall:

- Analytic: about 140 mm/s at t = 0.08 s
- Mesh: about 82 mm/s at t = 0.08 s

So the mesh rebound remains about 41% lower than the analytical result.

Also, your attached result shows avg_contacts = 1 during impact, while my mesh test shows avg_contacts = 6. Is this the number of sphere-triangle contact pairs? 

Could the different contact count, mesh triangulation, or impact position explain the remaining rebound difference?

Could you share the exact minimal test setup used for your attached result, especially the mesh, initial sphere position, and MM_TO_M = 1e-1 change?

Thank you.

2026년 8월 13일 목요일 오후 4시 11분 7초 UTC+9에 Ruochun Zhang님이 작성:

[Ruochun Zhang]

That's a good catch. Actually, if you change the sphere's initial position to 0.013, 0.020, 0.012 (or 1000 times that, if in mm) in your script, I think you will see the results I showed. I forgot to mention that was what I used.

The difference is that, in your initial position, the sphere hits the mesh at a point where some 6 triangles share a vertex. DEME2 estimates each sphere--triangle pair's effect separately, so the added effect may be close but not the same as hitting a perfect plane. The initial position I gave ensures the sphere hits a triangle's inner part, then it behaves the same as a plane. DEME2 has no way around this limitation but in general, this should not change the more important macro behavior of reasonably sized granular systems.

DEME3 actually improves it by combining the effect of all triangles involved, so it should provide a better estimation of the effect of meshed walls, but it still won't be the same as analytical. You might want to consider this if it is super critical for you.

Thank you,

Ruochun

[Ruochun Zhang]

I invite you to test your example using deme3. It can be found here. I would suggest that you create a new conda environment and then pip install deme3. In my testing, the sphere--mesh contact is much more accurately estimated using deme3.

Note that deme3 is still being tested and this deme3 PyPI package can be seen as a preview version. The official releases are still under deme.

Thank you,

Ruochun

[박선호]

Using DEME3, the mesh-particle rebound issue has been substantially improved, so I created the model shown in the figure.

Particles are initially placed in the upper chamber. They accumulate there, pass through the small outlet/funnel, and fall into the lower collection box, where they should form a packed bed.

However, after particles reach the lower box, nonphysical behavior appears. Some particles that are already rolling or resting on the bottom suddenly accelerate upward and scatter through the chamber. 

This behavior remains even when the coefficient of restitution is set to zero, so it is likely caused by unstable sphere-triangle contact impulses near the lower wall geometry rather than physical rebound.

2026년 8월 16일 일요일 오전 3시 7분 52초 UTC+9에 Ruochun Zhang님이 작성:

[Ruochun Zhang]

Again, could you share this script and the related mesh files? I may have seen something similar b before, but not reliably and perhaps with the Python version only. So I don't think it's a force model issue, and it may be related to contact detection or Python porting. Your script should be very helpful in verifying it.

Thank you,

Ruochun

[박선호]

 Please check with attached files

2026년 8월 24일 월요일 오후 3시 38분 37초 UTC+9에 Ruochun Zhang님이 작성:

[Ruochun Zhang]

I think this is an interesting observation, and the cause is multifaceted.

The divergence in DEME3 appears to be caused by merging sphere–triangle contacts across large angles on the mesh surface. Because the box has inward-facing normals, its corners are effectively highly concave regions. When a sphere hits a corner, contacting two or three mutually perpendicular walls, the solver attempts to combine the effects of those facets. This can produce a nonphysical effective contact point and destabilize the simulation.

This also explains why your script works in DEME2: DEME2 treats each triangle as an independent patch and does not merge the effects of adjacent facets.

To fix this in DEME3, call:

mesh.SplitIntoConvexPatches(45.0)

after loading each mesh. This separates triangles across surface angles greater than 45°, preventing contacts from being merged across sharp corners. You can think of it as dividing the mesh into smaller, better-behaved contact patches. This should resolve the problem.

Alternatively, call:

solver.SetDEME2MeshBehavior()

before loading any meshes. This restores the DEME2 behavior of treating every triangle as an individual patch. However, you might as well use DEME2 if you do this.

Both fixes require updating DEME3 to version 3.0.8. We are also investigating whether there is a more elegant, automatic and general solution for contacts in highly concave regions.

As a side note, I noticed that your script changes gravity after initialization. In older versions of DEME, including DEME2 and DEME3 3.0.7 or earlier, changing gravity alone does not immediately update the initialized workers. You must call:

solver.UpdateSimParams()

after each change for it to take effect. Starting with DEME3 3.0.8, gravitational-acceleration changes are propagated immediately and take effect on the next simulation step.

Thank you,

Ruochun