Skip to content

Add 3D beam elasticity simulation #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
th-skam merged 4 commits into from
Apr 3, 2026
Merged

Add 3D beam elasticity simulation #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5dec331
Add 3D beam elasticity: SOFA, FreeFEM and mesh file
Fimache Apr 2, 2026
0b123a9
Reworked script to include both plane stress and strain comparisons
th-skam Apr 2, 2026
4862f2c
Add a gmsh script for 3D mesh generation of a beam
th-skam Apr 3, 2026
c7aeed0
fix gmsh code and apply comparison script
th-skam Apr 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
199 changes: 102 additions & 97 deletions examples/Freefem/2D-beam/compare_sofa_freefem.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,106 +54,111 @@ def _write_freefem_results(path, x, y, u_x, u_y):
filename=mesh_filename
)

# --- Run FreeFEM ---
# gmsh.initialize() corrupts PATH, breaking pyfreefem's stdbuf call
os.environ['PATH'] = path_before_gmsh
runner = FreeFemRunner("freefem_beam2d.edp")
results = runner.execute({
'youngModulus': young_modulus,
'poissonRatio': poisson_ratio,
'rhoMat': rho,
'gravity': gravity,
'meshFile': os.path.abspath(msh_path).replace('\\', '/'),
})
x_ff = results['xcoords']
y_ff = results['ycoords']
dofs = results['uu[]'] # interleaved [uu_0, vv_0, uu_1, vv_1, ...]
u_x_ff = dofs[0::2]
u_y_ff = dofs[1::2]

os.makedirs("results", exist_ok=True)
_write_freefem_results(os.path.join("results", "freefem_results.txt"), x_ff, y_ff, u_x_ff, u_y_ff)

## --- Run SOFA ---
pos_initial, u_sofa = sofaRun(
height=height,
young_modulus=young_modulus,
poisson_ratio=poisson_ratio,
rho=rho,
gravity=gravity,
mesh_file=os.path.abspath(msh_path).replace('\\', '/')
)
x_sofa = pos_initial[:, 0]
y_sofa = pos_initial[:, 1]
u_x_sofa = u_sofa[:, 0]
u_y_sofa = u_sofa[:, 1]

## --- Compare Results ---
err_ux = u_x_sofa - u_x_ff
err_uy = u_y_sofa - u_y_ff

with open("results/comparison_results.txt", 'w') as f:
header = (f"{'x':>12} {'y':>12} {'ux_ff':>12} {'ux_sofa':>12}"
f" {'uy_ff':>12} {'uy_sofa':>12} {'err_ux':>12} {'err_uy':>12}")
f.write(header + "\n")
f.write("-" * len(header) + "\n")
for vals in zip(x_ff, y_ff, u_x_ff, u_x_sofa, u_y_ff, u_y_sofa, err_ux, err_uy):
f.write(" ".join(f"{v:12.6f}" for v in vals) + "\n")
f.write("\n")
f.write("L2 error norms\n")
f.write("-" * 40 + "\n")
f.write(f" ||ux_sofa - ux_ff||_2 = {_l2(u_x_sofa, u_x_ff):.6e}\n")
f.write(f" ||uy_sofa - uy_ff||_2 = {_l2(u_y_sofa, u_y_ff):.6e}\n")

## --- Plots ---
triang = mtri.Triangulation(x_ff, y_ff)

# Shared colour ranges so SOFA and FreeFEM panels are directly comparable
vmin_ux = min(u_x_ff.min(), u_x_sofa.min())
vmax_ux = max(u_x_ff.max(), u_x_sofa.max())
vmin_uy = min(u_y_ff.min(), u_y_sofa.min())
vmax_uy = max(u_y_ff.max(), u_y_sofa.max())

# Fig 1 — displacement field: FreeFEM (top row) vs SOFA (bottom row)
fig, axes = plt.subplots(2, 2, figsize=(14, 6), constrained_layout=True)
fig.suptitle("Displacement field — FreeFEM (top) vs SOFA (bottom)")
for row, (label, ux, uy) in enumerate([
("FreeFEM", u_x_ff, u_y_ff),
("SOFA", u_x_sofa, u_y_sofa),
]):
tcf = axes[row, 0].tricontourf(triang, ux, levels=20, vmin=vmin_ux, vmax=vmax_ux)
axes[row, 0].set_title(f"{label} — $u_x$")
axes[row, 0].set_aspect('equal')
fig.colorbar(tcf, ax=axes[row, 0])

tcf = axes[row, 1].tricontourf(triang, uy, levels=20, vmin=vmin_uy, vmax=vmax_uy)
axes[row, 1].set_title(f"{label} — $u_y$")
axes[row, 1].set_aspect('equal')
fig.colorbar(tcf, ax=axes[row, 1])
fig.savefig("results/plot_displacement_field.png", dpi=150)
plt.close(fig)

# Fig 2 — deformed beams overlaid
fig, ax = plt.subplots(figsize=(14, 4))
ax.scatter(x_ff + u_x_ff, y_ff + u_y_ff, s=4, label="FreeFEM", alpha=0.7)
ax.scatter(x_sofa + u_x_sofa, y_sofa + u_y_sofa, s=4, label="SOFA", alpha=0.7, marker='x')
ax.set_aspect('equal')
ax.set_title("Deformed beam — FreeFEM vs SOFA")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend()
fig.savefig("results/plot_deformed_beams.png", dpi=150)
plt.close(fig)

# Figs 3 & 4 — nodal error in x and y
for comp, err in [("ux", err_ux), ("uy", err_uy)]:
for formulation2d in ["planeStrain", "planeStress"]:
print(f"Running comparison simulations for {formulation2d}")
# --- Run FreeFEM ---
# gmsh.initialize() corrupts PATH, breaking pyfreefem's stdbuf call
os.environ['PATH'] = path_before_gmsh
runner = FreeFemRunner("freefem_beam2d.edp")
results = runner.execute({
'youngModulus': young_modulus,
'poissonRatio': poisson_ratio,
'rhoMat': rho,
'gravity': gravity,
'meshFile': os.path.abspath(msh_path).replace('\\', '/'),
'formulation' : formulation2d
})
x_ff = results['xcoords']
y_ff = results['ycoords']
dofs = results['uu[]'] # interleaved [uu_0, vv_0, uu_1, vv_1, ...]
u_x_ff = dofs[0::2]
u_y_ff = dofs[1::2]

os.makedirs("results", exist_ok=True)
_write_freefem_results(os.path.join("results", f"freefem_{formulation2d}_results.txt")
, x_ff, y_ff, u_x_ff, u_y_ff)

## --- Run SOFA ---
pos_initial, u_sofa = sofaRun(
height=height,
young_modulus=young_modulus,
poisson_ratio=poisson_ratio,
rho=rho,
gravity=gravity,
mesh_file=os.path.abspath(msh_path).replace('\\', '/'),
formulation=formulation2d
)
x_sofa = pos_initial[:, 0]
y_sofa = pos_initial[:, 1]
u_x_sofa = u_sofa[:, 0]
u_y_sofa = u_sofa[:, 1]

## --- Compare Results ---
err_ux = u_x_sofa - u_x_ff
err_uy = u_y_sofa - u_y_ff

with open(f"results/comparison_{formulation2d}_results.txt", 'w') as f:
header = (f"{'x':>12} {'y':>12} {'ux_ff':>12} {'ux_sofa':>12}"
f" {'uy_ff':>12} {'uy_sofa':>12} {'err_ux':>12} {'err_uy':>12}")
f.write(header + "\n")
f.write("-" * len(header) + "\n")
for vals in zip(x_ff, y_ff, u_x_ff, u_x_sofa, u_y_ff, u_y_sofa, err_ux, err_uy):
f.write(" ".join(f"{v:12.6f}" for v in vals) + "\n")
f.write("\n")
f.write("L2 error norms\n")
f.write("-" * 40 + "\n")
f.write(f" ||ux_sofa - ux_ff||_2 = {_l2(u_x_sofa, u_x_ff):.6e}\n")
f.write(f" ||uy_sofa - uy_ff||_2 = {_l2(u_y_sofa, u_y_ff):.6e}\n")

## --- Plots ---
triang = mtri.Triangulation(x_ff, y_ff)

# Shared colour ranges so SOFA and FreeFEM panels are directly comparable
vmin_ux = min(u_x_ff.min(), u_x_sofa.min())
vmax_ux = max(u_x_ff.max(), u_x_sofa.max())
vmin_uy = min(u_y_ff.min(), u_y_sofa.min())
vmax_uy = max(u_y_ff.max(), u_y_sofa.max())

# Fig 1 — displacement field: FreeFEM (top row) vs SOFA (bottom row)
fig, axes = plt.subplots(2, 2, figsize=(14, 6), constrained_layout=True)
fig.suptitle("Displacement field — FreeFEM (top) vs SOFA (bottom)")
for row, (label, ux, uy) in enumerate([
("FreeFEM", u_x_ff, u_y_ff),
("SOFA", u_x_sofa, u_y_sofa),
]):
tcf = axes[row, 0].tricontourf(triang, ux, levels=20, vmin=vmin_ux, vmax=vmax_ux)
axes[row, 0].set_title(f"{label} — $u_x$")
axes[row, 0].set_aspect('equal')
fig.colorbar(tcf, ax=axes[row, 0])

tcf = axes[row, 1].tricontourf(triang, uy, levels=20, vmin=vmin_uy, vmax=vmax_uy)
axes[row, 1].set_title(f"{label} — $u_y$")
axes[row, 1].set_aspect('equal')
fig.colorbar(tcf, ax=axes[row, 1])
fig.savefig(f"results/plot_{formulation2d}_displacement_field.png", dpi=150)
plt.close(fig)

# Fig 2 — deformed beams overlaid
fig, ax = plt.subplots(figsize=(14, 4))
tcf = ax.tricontourf(triang, err, levels=20, cmap='RdBu_r')
fig.colorbar(tcf, ax=ax, label=f"$u_{comp[1]}$ error (SOFA − FreeFEM)")
ax.set_title(f"Nodal error — $u_{comp[1]}$")
ax.scatter(x_ff + u_x_ff, y_ff + u_y_ff, s=4, label="FreeFEM", alpha=0.7)
ax.scatter(x_sofa + u_x_sofa, y_sofa + u_y_sofa, s=4, label="SOFA", alpha=0.7, marker='x')
ax.set_aspect('equal')
ax.set_title("Deformed beam — FreeFEM vs SOFA")
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.savefig(f"results/plot_error_{comp}.png", dpi=150)
ax.legend()
fig.savefig(f"results/plot_{formulation2d}_deformed_beams.png", dpi=150)
plt.close(fig)

# Figs 3 & 4 — nodal error in x and y
for comp, err in [("ux", err_ux), ("uy", err_uy)]:
fig, ax = plt.subplots(figsize=(14, 4))
tcf = ax.tricontourf(triang, err, levels=20, cmap='RdBu_r')
fig.colorbar(tcf, ax=ax, label=f"$u_{comp[1]}$ error (SOFA − FreeFEM)")
ax.set_title(f"Nodal error — $u_{comp[1]}$")
ax.set_aspect('equal')
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.savefig(f"results/plot_error_{formulation2d}_{comp}.png", dpi=150)
plt.close(fig)

11 changes: 8 additions & 3 deletions examples/Freefem/2D-beam/freefem_beam2d.edp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,10 @@ DEFAULT (rhoMat, "7.85e-6")
DEFAULT (gravity, "9810")
DEFAULT (poissonRatio, "0.3")
DEFAULT (meshFile, "")
DEFAULT (formulation, "planeStrain")

real E = $youngModulus;
real sigma = $poissonRatio;
real nu = $poissonRatio;
real rhoMat = $rhoMat;
real gravity = $gravity;
real F = -rhoMat * gravity;
Expand All @@ -22,8 +23,12 @@ real sqrt2 = sqrt(2.);
macro epsilon(u1,u2) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u,v) (dx(u)+dy(v)) // EOM

real mu = E / (2*(1 + sigma));
real lambda = E*sigma / ((1 + sigma)*(1 - 2*sigma));
real mu = E / (2*(1 + nu));
IFEQ (formulation, "planeStrain")
real lambda = E*nu / ((1 + nu)*(1 - 2*nu)); // 3D but z is constrained
ELSE
real lambda = E*nu / (1 - nu*nu); // Pure 2D
ENDIF

solve Elasticity([uu, vv], [w, s], solver=sparsesolver)
= int2d(th)(
Expand Down
33 changes: 23 additions & 10 deletions examples/Freefem/2D-beam/sofa_beam2d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,16 @@ def onAnimateEndEvent(self, event):
f.write(f"{x:12.6f} {y:12.6f} {ux:12.6f} {uy:12.6f}\n")


def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio, rho, gravity):
def create_scene_args(rootNode
, mesh_file
, height
, young_modulus
, poisson_ratio
, rho
, gravity
, formulation):
dim_template = "Vec3d" if formulation == "planeStrain" else "Vec2d"

requiredPlugins = [
"Elasticity",
"Sofa.Component.Constraint.Projective",
Expand Down Expand Up @@ -85,7 +94,7 @@ def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio,
, filename=mesh_file)
dofs = Beam.addObject('MechanicalObject'
, name="dofs"
, template="Vec3d"
, template=dim_template
, src="@loader")

with Beam.addChild('triangles') as Triangles:
Expand All @@ -94,15 +103,16 @@ def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio,
, src="@../loader")
Triangles.addObject('TriangleSetTopologyModifier')
Triangles.addObject('TriangleSetGeometryAlgorithms'
, template="Vec3d"
, template=dim_template
, drawTriangles=True)
Triangles.addObject('MeshMatrixMass'
, name="mass"
, template=f"{dim_template},{dim_template}"
, massDensity=rho
, topology="@topology")
Triangles.addObject('LinearSmallStrainFEMForceField'
, name="FEM"
, template="Vec3d"
, template=dim_template
, youngModulus=young_modulus
, poissonRatio=poisson_ratio
, topology="@topology")
Expand All @@ -111,18 +121,18 @@ def create_scene_args(rootNode, mesh_file, height, young_modulus, poisson_ratio,
eps = 1e-5
Beam.addObject('BoxROI'
, name="fixed_roi"
, template="Vec3d"
, template=dim_template
, box=[-eps, -eps, -1.0, eps, height + eps, 1.0]
, drawBoxes=True)
Beam.addObject('FixedProjectiveConstraint'
, template="Vec3d"
, template=dim_template
, indices="@fixed_roi.indices")

os.makedirs(RESULTS_DIR, exist_ok=True)
exporter = rootNode.addObject(
DisplacementExporter(
dofs_node = dofs,
output_file = os.path.join(RESULTS_DIR, "sofa_results.txt"),
output_file = os.path.join(RESULTS_DIR, f"sofa_{formulation}_results.txt"),
name = "exportCtrl"
)
)
Expand All @@ -140,19 +150,21 @@ def createScene(rootNode):
, rho = float(cfg["rho"])
, gravity = float(cfg["gravity"])
, mesh_file = os.path.join(RESULTS_DIR, cfg.get("meshfile", "beam2d.msh"))
, formulation = "planeStrain"
)
return rootNode


def sofaRun(height, young_modulus, poisson_ratio, rho, gravity, mesh_file=None):
def sofaRun(height, young_modulus, poisson_ratio, rho, gravity, mesh_file, formulation):
root = Sofa.Core.Node("root")
_, exporter = create_scene_args(root
, mesh_file = mesh_file
, height = height
, height = height
, young_modulus= young_modulus
, poisson_ratio= poisson_ratio
, rho = rho
, gravity = gravity)
, gravity = gravity
, formulation = formulation)
Sofa.Simulation.init(root)
Sofa.Simulation.animate(root, root.dt.value)
return exporter.pos_initial, exporter.u
Expand All @@ -169,4 +181,5 @@ def sofaRun(height, young_modulus, poisson_ratio, rho, gravity, mesh_file=None):
, rho = float(cfg["rho"])
, gravity = float(cfg["gravity"])
, mesh_file = os.path.join(RESULTS_DIR, cfg.get("meshfile", "beam2d.msh"))
, formulation = "planeStrain"
)
Loading
Loading