import gmsh
class VerticalSphereStack:
def __init__(self, number, radius, spacing, box):
self.number = number
self.radius = radius
self.spacing = spacing
self.box = box # [x_min, y_min, z_min, x_max, y_max, z_max]
self.objects = []
def generate_positions(self):
y_start = self.box[1] + self.radius # Start from bottom boundary with radius offset
for i in range(self.number):
y_position = y_start + i * (2 * self.radius - self.spacing)
self.objects.append({"x": 0, "y": y_position, "z": 0, "radius": self.radius})
def generate_geometry(self):
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("Vertical Sphere Stack")
# Add spheres
for i, obj in enumerate(self.objects):
gmsh.model.occ.addSphere(obj["x"], obj["y"], obj["z"], obj["radius"], tag=i + 1)
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.02)
try:
gmsh.model.mesh.generate(3)
except Exception as e:
print(f"Error during mesh generation: {e}")
gmsh.write("finer.msh")
gmsh.write("vertical_sphere_stack.step")
gmsh.finalize()
def run(self):
print("Generating sphere positions...")
self.generate_positions()
print("Generating Gmsh geometry...")
self.generate_geometry()
print("Done! Files saved as 'vertical_sphere_stack.msh' and 'vertical_sphere_stack.step'.")
# Parameters
number_of_spheres = 2
particle_radius = 0.1
spacing_between_spheres = 0.001 # Small gap between spheres
bounding_box = [-1, -1, -1, 1, 1, 1] # [x_min, y_min, z_min, x_max, y_max, z_max]
sphere_stack = VerticalSphereStack(number_of_spheres, particle_radius, spacing_between_spheres, bounding_box)
sphere_stack.run()
Minimal example for compressing particle
from dolfin import *
import numpy as np
import meshio
# Read mesh
mesh_from_file = meshio.read("finer.msh")
# Extract tetrahedral elements
tetra_cells = mesh_from_file.cells_dict.get("tetra", None)
if tetra_cells is None or len(tetra_cells) == 0:
raise ValueError("No tetrahedral elements found in the mesh.")
# Write tetrahedral mesh to XDMF
meshio.write_points_cells(
"finer.xdmf",
points=mesh_from_file.points,
cells=[("tetra", tetra_cells)]
)
print("Mesh successfully written to finer.xdmf")
# Load the Mesh into FEniCS
mesh = Mesh()
with XDMFFile("finer.xdmf") as infile:
infile.read(mesh)
V = VectorFunctionSpace(mesh, "P", 1)
# Extract region markers
if "gmsh:physical" not in mesh_from_file.cell_data:
raise ValueError("Mesh file does not contain 'gmsh:physical'. Check if the mesh was generated correctly.")
physical_region_data = mesh_from_file.cell_data["gmsh:physical"][0]
# **Check and print available region IDs**
unique_regions = np.unique(physical_region_data)
print(f"Detected physical region IDs in the mesh: {unique_regions}")
# Define boundary facets based on region markers
facet_regions = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
# Assign facets to boundary regions
for i, cell in enumerate(mesh_from_file.cells[-1].data): # Extract boundary facets
region_id = physical_region_data[i]
facet_regions.array()[i] = region_id # Assign detected region ID
# Save boundary markers for debugging
File("boundary_markers.pvd") << facet_regions
# **Identify Region1 (Fixed) and Region2 (Compression)**
region1_id = 13 if 13 in unique_regions else min(unique_regions) # Fully fixed
region2_id = 14 if 14 in unique_regions else max(unique_regions) # Compressive displacement
print(f"Applying BCs on Region1 ID: {region1_id} (Fully Fixed), Region2 ID: {region2_id} (Compression towards Region1)")
# **Ensure Region1 is Fully Fixed**
# Directly use facet_regions to detect Region1 boundary facets
bc_region1 = DirichletBC(V, Constant((0.0, 0.0, 0.0)), facet_regions, region1_id)
# **Apply compressive displacement on Region2 (Moving towards Region1)**
alpha = 0.01 # Compression rate
t_start, t_end, dt = 0.0, 100.0, 1
time_steps = np.arange(t_start, t_end + dt, dt)
bc_region2 = DirichletBC(V, Expression(("0.0", "-alpha * t", "0.0"), alpha=alpha, t=t_start, degree=1), facet_regions, region2_id)
# Apply boundary conditions
bcs = [bc_region1, bc_region2]
# Time-stepping loop
u_trial, v = TrialFunction(V), TestFunction(V)
u_solution = Function(V)
grouped_displacement_file = File("outputs/grouped_displacement.pvd")
# Material properties
E, nu = 1.0, 0.3
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# Define elasticity formulation
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lmbda * div(u) * Identity(3) + 2 * mu * epsilon(u)
# **Contact Mechanics Using Hertzian Model**
R = 0.2 # Contact radius
gap_tol = 1e-3 # Allowable gap before contact force applies
# Define gap function (how much Region2 moves towards Region1)
gap = Expression("x[1] < gap_tol ? gap_tol - x[1] : 0", gap_tol=gap_tol, degree=1)
# Define Hertzian contact force: Only applies if the gap is negative (contact occurs)
def hertzian_force(delta):
return conditional(gt(delta, 0), (4/3) * E * np.sqrt(R) * delta**(3/2), 0)
contact_force = hertzian_force(gap)
# **Final Weak Form of the Problem**
a = inner(sigma(u_trial), epsilon(v)) * dx
L = dot(Constant((0, 0, 0)), v) * dx + contact_force * v[1] * dx # Contact force added
for t in time_steps:
bc_region2.set_value(Expression(("0.0", "-alpha * t", "0.0"), alpha=alpha, t=t, degree=1))
print(f"Solving for time t = {t}")
solve(a == L, u_solution, bcs)
# **Ensure Proper Visualization**
u_projected = project(u_solution, V) # Ensure displacement field is fully computed
u_projected.rename("Displacement", "u") # Rename to avoid ParaView issues
grouped_displacement_file << (u_projected, t) # Save output for visualization
print("Compression simulation with contact completed.")
Note: Regions i defined in Gmsh itself, so i have used regions in code.