I'm new to the field, so I apologize if this is a basic question. I'm working with PyDrake's DeformableModel to simulate a soft joint between a rigid column and a rigid leg (both URDF files). The joint is modeled as a deformable cylinder through a volumetric thetraedrical mesh (.vtk), but when I run the simulation, it has no degrees of freedom (DOF). It has dimensions (0.06m,0.06m,0.08m). This is the first issue.

I try to attach the joint over the colum by using deformable_model.AddFixedConstraint() function, however even before constraint it has no degrees of freedom. Plus, when i constrain it choosing a a box as a shape no constraint is added to the model. The joint must place inside the upper face of the column (also a cylinder without holes) for 0.02 m, which i'm doing with the constraint.

I'll share the code for clarity:

import numpy as np
import os
import pydrake

from pydrake.all import *


# Check where we are
#print("Current directory:", os.getcwd())

#fixing timing 
simulation_time = 8.0   # Desired duration of the simulation [s].
realtime_rate = 1.0     # Desired real time rate.
time_step = 1e-2   

     # Discrete time step for the system [s]. Must be positive.
render_flag = True

## ---------------------------------------------------------------------------- COORDINATE ISSUE ----------------------------------------------------------------------------------
#   1. Rigid Trasform takes CENTER OF MASS coordinate, so work with them only
#   2. The joint is placed for 0.01 m inside the leg and for the same quantity inside the column
#   3. All of the coordinate I give here are referred to the external system of reference, AKA the ground. Everything is considering the penetrating joint and in meters.

# column
column_height = 0.20
z_column_center = 0.1

# joint 
joint_height = 0.08
joint_z_center = -0.02 #already at 0.24 in the vtk files
joint_radius = 0.055

# leg
leg_height = 0.675
leg_z_center = 0.5975

## --------------------------------------------------------------------------------------- MAIN PLANT ----------------------------------------------------------------------------

# Initiliazing the model
builder = DiagramBuilder()           #creating plant

plant_config = MultibodyPlantConfig()
plant_config.time_step = time_step

plant, scene_graph = AddMultibodyPlant(plant_config, builder)

#Add a renderer to render the inside dot pattern of the bubble gripper.
#Currently (April 2024), deformable rendering is only supported by
#RenderEngineGl.
if render_flag:
    scene_graph.AddRenderer("gl_renderer", MakeRenderEngineGl(RenderEngineGlParams()))

#Minimum required proximity properties for rigid bodies to interact with
#deformable bodies.
#1. A valid Coulomb friction coefficient, and
#2. A resolution hint. (Rigid bodies need to be tessellated so that collision
#queries can be performed against deformable geometries.) 


rigid_proximity_props = ProximityProperties()
surface_friction = CoulombFriction(1.15, 1.15)
resolution_hint = 0.01;
AddContactMaterial(friction = surface_friction,  properties = rigid_proximity_props)
rigid_proximity_props.AddProperty("hydroelastic", "resolution_hint", resolution_hint)


#Rigid parts
parser = Parser(plant)              # Parser for URDF files

# Adjust paths to correctly point to URDF files
base_column_path = "colored_oneleg_column.urdf"
leg_path = "colored_oneleg_leg.urdf"

# Verify file paths before loading URDF files
if not os.path.exists(base_column_path):
    raise FileNotFoundError(f"File not found: {base_column_path}")

if not os.path.exists(leg_path):
    raise FileNotFoundError(f"File not found: {leg_path}")


column_models = parser.AddModels(base_column_path)[0]

# Get the column instance and its frame
colonna_model_instance = column_models  # Ensure you're using the correct instance
body_base = plant.GetBodyByName("Link_colonna", colonna_model_instance)
column_frame = body_base.body_frame()

# Weld the column to the world frame
plant.WeldFrames(
    plant.world_frame(), 
    column_frame, 
    RigidTransform([0, 0, z_column_center])
)



## ---------------------------------------------------------------------------- DEFORMABLE PART (JOINT) ---------------------------------------------------------------------------------------

# Material properties for the elastic joint (steel)
density_ms = 7800          # kg/m³
youngs_modulus_ms = 200e9  # Young's Modulus for Maraging steel (Pa)
poisson_ratio_ms = 0.3     # Poisson's ratio for steel
beta = 0.01                # Stiffness damping coefficient for the deformable body [1/s].

# Assign the  material model with linear corotated material model
#documentation: .multibody.fem.html?highlight=deformablebodyconfig#pydrake.multibody.fem.MaterialModel
material_model_joint = MaterialModel.kLinearCorotated

#Geometry
joint_path = "volumetric_vtk_joint_no_triangles.vtk"  # Update with the actual mesh file path
if not os.path.exists(joint_path):
    raise FileNotFoundError(f"Mesh file not found: {joint_path}")

#pose the joint in world frame so that it is place over the column for 0.01
X_WB = RigidTransform(RotationMatrix(), [0, 0, joint_z_center])

# geometry instance for the joint: define mesh and place it
joint_geometry = GeometryInstance(X_WB,  Mesh("volumetric_vtk_joint_no_triangles.vtk"),  "deformable_joint")

#jont properties
joint_deformable_model = DeformableModel(plant)
config = DeformableBodyConfig()
config.set_youngs_modulus(youngs_modulus=youngs_modulus_ms)
config.set_poissons_ratio(poissons_ratio=poisson_ratio_ms)
config.set_stiffness_damping_coefficient(beta)
config.set_mass_density(mass_density=density_ms)
config.set_material_model(material_model_joint)


# Minimumly required proximity properties for deformable bodies:
# A valid Coulomb friction coefficient.
deformable_proximity_props = ProximityProperties()
AddContactMaterial(friction=surface_friction, properties=deformable_proximity_props)
joint_geometry.set_proximity_properties(deformable_proximity_props)

fem_resolution_hint = 0.3  # Controls FEM discretization
joint_deformable_model = plant.mutable_deformable_model()

#register the deformbale model
deformable_body_id = joint_deformable_model.RegisterDeformableBody(
    geometry_instance=joint_geometry,
    config=config,
    resolution_hint=fem_resolution_hint
)
deformable_model = joint_deformable_model
print(deformable_model.GetReferencePositions(deformable_body_id))

#Now we attach the bubble to the WSG finger using a fixed constraint. To do
#that, we specify a box geometry and put all vertices of the bubble geometry
#under fixed constraint with the rigid finger if they fall inside the box.
#Refer to DeformableModel::AddFixedConstraint for details. 


#joint in the column frame
X_JC = RigidTransform(RotationMatrix(), np.array([0,0, 0.08]))
#All vertices of the deformable bubble mesh inside this box will be subject 
box_column = Box(0.06, 0.06, 0.02)  

#attach to column (see attched sheet for the reasoning)
#deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) ) 

plant.Finalize()    #finalize plant

# Check the number of DOF and velocities for the deformable body
num_dof = plant.num_velocities()
num_joints = plant.num_joints()
num_frames = plant.num_frames()
num_pos = plant.num_positions()
 
print(f"Number of velocities: {num_dof}")
print(f"Number of joints: {num_joints}")
print(f"Number of frames: {num_frames}")
print(f"Number of positions: {num_pos}")
num_deformable_bodies = deformable_model.num_bodies()

print(f"Number of deformable bodies: {num_deformable_bodies}")



### ----------------------------------------------------------------------- SIMULATION AND RENDERING ------------------------------------------------------------
#builder
DrakeVisualizer.AddToBuilder(builder, scene_graph)
diagram = builder.Build()
diagram_context = diagram.CreateDefaultContext()

#simulator
simulator = Simulator(diagram, diagram_context)
mutable_root_context = simulator.get_mutable_context()
plant_context =diagram.GetMutableSubsystemContext(plant, mutable_root_context)
simulator.Initialize()

#timing
simulator.set_target_realtime_rate(realtime_rate)
simulator.AdvanceTo(simulation_time)

Output commenting the deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) ) function: [ 0.0005365 0.00544716 0.24999998 ... -0.00211268 -0.00385886 0.25932274] Number of velocities: 0 Number of joints: 1 Number of frames: 2 Number of positions: 0 Number of deformable bodies: 1

Output without commenting it: `[ 0.0005365 0.00544716 0.24999998 ... -0.00211268 -0.00385886 0.25932274] Traceback (most recent call last): File "/home/gaia/Scrivania/deformable_drake/deformable_joint.py", line 160, in deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) ) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ RuntimeError: No constraint has been added between deformable body with id 32 and rigid body with name Link_colonna. Remove the call to AddFixedConstraint() if this is intended.

`

I'm pretty sure it is not a coordinate problem, beacuse I've double checked by visualizing the model. I've also check the mesh and it is made of Tetrahedrons only. How can i debug it? Is there a guide for deformable models?

I'm running the code on Ubuntu 24.04.1

I'm new to the field, so I apologize if this is a basic question. I'm working with PyDrake's DeformableModel to simulate a soft joint between a rigid column and a rigid leg (both URDF files). The joint is modeled as a deformable cylinder through a volumetric thetraedrical mesh (.vtk), but when I run the simulation, it has no degrees of freedom (DOF). It has dimensions (0.06m,0.06m,0.08m). This is the first issue.

I try to attach the joint over the colum by using deformable_model.AddFixedConstraint() function, however even before constraint it has no degrees of freedom. Plus, when i constrain it choosing a a box as a shape no constraint is added to the model. The joint must place inside the upper face of the column (also a cylinder without holes) for 0.02 m, which i'm doing with the constraint.

I'll share the code for clarity:

import numpy as np
import os
import pydrake

from pydrake.all import *


# Check where we are
#print("Current directory:", os.getcwd())

#fixing timing 
simulation_time = 8.0   # Desired duration of the simulation [s].
realtime_rate = 1.0     # Desired real time rate.
time_step = 1e-2   

     # Discrete time step for the system [s]. Must be positive.
render_flag = True

## ---------------------------------------------------------------------------- COORDINATE ISSUE ----------------------------------------------------------------------------------
#   1. Rigid Trasform takes CENTER OF MASS coordinate, so work with them only
#   2. The joint is placed for 0.01 m inside the leg and for the same quantity inside the column
#   3. All of the coordinate I give here are referred to the external system of reference, AKA the ground. Everything is considering the penetrating joint and in meters.

# column
column_height = 0.20
z_column_center = 0.1

# joint 
joint_height = 0.08
joint_z_center = -0.02 #already at 0.24 in the vtk files
joint_radius = 0.055

# leg
leg_height = 0.675
leg_z_center = 0.5975

## --------------------------------------------------------------------------------------- MAIN PLANT ----------------------------------------------------------------------------

# Initiliazing the model
builder = DiagramBuilder()           #creating plant

plant_config = MultibodyPlantConfig()
plant_config.time_step = time_step

plant, scene_graph = AddMultibodyPlant(plant_config, builder)

#Add a renderer to render the inside dot pattern of the bubble gripper.
#Currently (April 2024), deformable rendering is only supported by
#RenderEngineGl.
if render_flag:
    scene_graph.AddRenderer("gl_renderer", MakeRenderEngineGl(RenderEngineGlParams()))

#Minimum required proximity properties for rigid bodies to interact with
#deformable bodies.
#1. A valid Coulomb friction coefficient, and
#2. A resolution hint. (Rigid bodies need to be tessellated so that collision
#queries can be performed against deformable geometries.) 


rigid_proximity_props = ProximityProperties()
surface_friction = CoulombFriction(1.15, 1.15)
resolution_hint = 0.01;
AddContactMaterial(friction = surface_friction,  properties = rigid_proximity_props)
rigid_proximity_props.AddProperty("hydroelastic", "resolution_hint", resolution_hint)


#Rigid parts
parser = Parser(plant)              # Parser for URDF files

# Adjust paths to correctly point to URDF files
base_column_path = "colored_oneleg_column.urdf"
leg_path = "colored_oneleg_leg.urdf"

# Verify file paths before loading URDF files
if not os.path.exists(base_column_path):
    raise FileNotFoundError(f"File not found: {base_column_path}")

if not os.path.exists(leg_path):
    raise FileNotFoundError(f"File not found: {leg_path}")


column_models = parser.AddModels(base_column_path)[0]

# Get the column instance and its frame
colonna_model_instance = column_models  # Ensure you're using the correct instance
body_base = plant.GetBodyByName("Link_colonna", colonna_model_instance)
column_frame = body_base.body_frame()

# Weld the column to the world frame
plant.WeldFrames(
    plant.world_frame(), 
    column_frame, 
    RigidTransform([0, 0, z_column_center])
)



## ---------------------------------------------------------------------------- DEFORMABLE PART (JOINT) ---------------------------------------------------------------------------------------

# Material properties for the elastic joint (steel)
density_ms = 7800          # kg/m³
youngs_modulus_ms = 200e9  # Young's Modulus for Maraging steel (Pa)
poisson_ratio_ms = 0.3     # Poisson's ratio for steel
beta = 0.01                # Stiffness damping coefficient for the deformable body [1/s].

# Assign the  material model with linear corotated material model
#documentation: https://drake.mit.edu/pydrake/pydrake.multibody.fem.html?highlight=deformablebodyconfig#pydrake.multibody.fem.MaterialModel
material_model_joint = MaterialModel.kLinearCorotated

#Geometry
joint_path = "volumetric_vtk_joint_no_triangles.vtk"  # Update with the actual mesh file path
if not os.path.exists(joint_path):
    raise FileNotFoundError(f"Mesh file not found: {joint_path}")

#pose the joint in world frame so that it is place over the column for 0.01
X_WB = RigidTransform(RotationMatrix(), [0, 0, joint_z_center])

# geometry instance for the joint: define mesh and place it
joint_geometry = GeometryInstance(X_WB,  Mesh("volumetric_vtk_joint_no_triangles.vtk"),  "deformable_joint")

#jont properties
joint_deformable_model = DeformableModel(plant)
config = DeformableBodyConfig()
config.set_youngs_modulus(youngs_modulus=youngs_modulus_ms)
config.set_poissons_ratio(poissons_ratio=poisson_ratio_ms)
config.set_stiffness_damping_coefficient(beta)
config.set_mass_density(mass_density=density_ms)
config.set_material_model(material_model_joint)


# Minimumly required proximity properties for deformable bodies:
# A valid Coulomb friction coefficient.
deformable_proximity_props = ProximityProperties()
AddContactMaterial(friction=surface_friction, properties=deformable_proximity_props)
joint_geometry.set_proximity_properties(deformable_proximity_props)

fem_resolution_hint = 0.3  # Controls FEM discretization
joint_deformable_model = plant.mutable_deformable_model()

#register the deformbale model
deformable_body_id = joint_deformable_model.RegisterDeformableBody(
    geometry_instance=joint_geometry,
    config=config,
    resolution_hint=fem_resolution_hint
)
deformable_model = joint_deformable_model
print(deformable_model.GetReferencePositions(deformable_body_id))

#Now we attach the bubble to the WSG finger using a fixed constraint. To do
#that, we specify a box geometry and put all vertices of the bubble geometry
#under fixed constraint with the rigid finger if they fall inside the box.
#Refer to DeformableModel::AddFixedConstraint for details. 


#joint in the column frame
X_JC = RigidTransform(RotationMatrix(), np.array([0,0, 0.08]))
#All vertices of the deformable bubble mesh inside this box will be subject 
box_column = Box(0.06, 0.06, 0.02)  

#attach to column (see attched sheet for the reasoning)
#deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) ) 

plant.Finalize()    #finalize plant

# Check the number of DOF and velocities for the deformable body
num_dof = plant.num_velocities()
num_joints = plant.num_joints()
num_frames = plant.num_frames()
num_pos = plant.num_positions()
 
print(f"Number of velocities: {num_dof}")
print(f"Number of joints: {num_joints}")
print(f"Number of frames: {num_frames}")
print(f"Number of positions: {num_pos}")
num_deformable_bodies = deformable_model.num_bodies()

print(f"Number of deformable bodies: {num_deformable_bodies}")



### ----------------------------------------------------------------------- SIMULATION AND RENDERING ------------------------------------------------------------
#builder
DrakeVisualizer.AddToBuilder(builder, scene_graph)
diagram = builder.Build()
diagram_context = diagram.CreateDefaultContext()

#simulator
simulator = Simulator(diagram, diagram_context)
mutable_root_context = simulator.get_mutable_context()
plant_context =diagram.GetMutableSubsystemContext(plant, mutable_root_context)
simulator.Initialize()

#timing
simulator.set_target_realtime_rate(realtime_rate)
simulator.AdvanceTo(simulation_time)

Output commenting the deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) ) function: [ 0.0005365 0.00544716 0.24999998 ... -0.00211268 -0.00385886 0.25932274] Number of velocities: 0 Number of joints: 1 Number of frames: 2 Number of positions: 0 Number of deformable bodies: 1

Output without commenting it: `[ 0.0005365 0.00544716 0.24999998 ... -0.00211268 -0.00385886 0.25932274] Traceback (most recent call last): File "/home/gaia/Scrivania/deformable_drake/deformable_joint.py", line 160, in deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) ) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ RuntimeError: No constraint has been added between deformable body with id 32 and rigid body with name Link_colonna. Remove the call to AddFixedConstraint() if this is intended.

`

I'm pretty sure it is not a coordinate problem, beacuse I've double checked by visualizing the model. I've also check the mesh and it is made of Tetrahedrons only. How can i debug it? Is there a guide for deformable models?

I'm running the code on Ubuntu 24.04.1

• drake
• finite-element-analysis

1 Answer 1

The fact that your print out shows "Number of deformable bodies: 1" indicates that the deformable body has been successfully added.

I'd recommend that you review this line:

deformable_model.AddFixedConstraint(deformable_body_id, body_base, X_JC, box_column, RigidTransform(RotationMatrix(), [0, 0, 0.085]) )

What you are trying to do here is to

1. put the deformable body in the body_base frame, with the pose X_JC, and,
2. put the box geometry box_column in the body_base frame, with pose RigidTransform([0, 0, 0.085]), (let me name that X_BG` where B stands for body base and G stands for geometry), and finally
3. Constrain all vertices of the deformable body that are inside the box_column to the body_base.

The error message indicates that if you prescribe poses X_JC and X_BG like that, NONE of the deformable body's vertices is falling inside the box geometry.

You should review the these poses.