Hello guys,
I am a beginner to fenicsx, and I am trying to wrap an elastic band around a circle, setting the boundary conditions for the left and right edges to meet at the top. I am either doing it wrong, or pyvista can’t visualize in the negative direction
This first code snippet generates the mesh with gmesh
import gmsh
import numpy as np
gmsh.initialize()
gmsh.model.add("forearm_electrode")
#Define forearm geometry (mm)
r = 45 # Radius of forearm
R_bone = r * np.sqrt(0.107)
R_muscle = r * np.sqrt(0.612 + 0.107)
R_fat = r * np.sqrt(0.223 + 0.612 + 0.107)
R_skin = r
# Band dimensions (mm)
L0 = 2 * np.pi * 30
t = 1 # Band thickness
Ed = 7
Et = 0.02
count = 16
IED = (L0 - (Ed * count))/count
bone = gmsh.model.occ.addDisk(0.0, 0.0, 0.0, R_bone, R_bone)
muscle = gmsh.model.occ.addDisk(0.0, 0.0, 0.0, R_muscle, R_muscle)
fat = gmsh.model.occ.addDisk(0.0, 0.0, 0.0, R_fat, R_fat)
skin = gmsh.model.occ.addDisk(0.0, 0.0, 0.0, R_skin, R_skin)
gmsh.model.occ.synchronize()
# Subtracting to get annuli
skin_ring = gmsh.model.occ.cut([(2, skin)], [(2, fat)], removeObject=True,
removeTool=False)
fat_ring = gmsh.model.occ.cut([(2, fat)], [(2, muscle)], removeObject=True,
removeTool=False)
muscle_ring = gmsh.model.occ.cut([(2, muscle)], [(2, bone)], removeObject=True,
removeTool=False)
gmsh.model.occ.synchronize()
print("Skin ring:", skin_ring)
# Collect resulting surface tags
skin_surfaces = [ent[1] for ent in skin_ring[0]]
fat_surfaces = [ent[1] for ent in fat_ring[0]]
muscle_surfaces = [ent[1] for ent in muscle_ring[0]]
bone_surfaces = [bone]
print("Skin surfaces:", skin_surfaces)
# Creating physical groups for surfaces
skin_tag = gmsh.model.addPhysicalGroup(2, skin_surfaces, name="Skin")
fat_tag = gmsh.model.addPhysicalGroup(2, fat_surfaces, name="Fat")
muscle_tag = gmsh.model.addPhysicalGroup(2, muscle_surfaces, name="Muscle")
bone_tag = gmsh.model.addPhysicalGroup(2, bone_surfaces, name="Bone")
# Creating physical groups for boundaries
bone_boundaries = gmsh.model.getBoundary([(2, tag) for tag in bone_surfaces], oriented=False)
muscle_boundaries = gmsh.model.getBoundary([(2, tag) for tag in muscle_surfaces], oriented=False)
fat_boundaries = gmsh.model.getBoundary([(2, tag) for tag in fat_surfaces], oriented=False)
skin_boundaries = gmsh.model.getBoundary([(2, tag) for tag in skin_surfaces], oriented=False)
gmsh.model.occ.synchronize()
print("Bone boundary curves:", bone_boundaries)
print("Muscle boundary curves:", muscle_boundaries)
print("Fat boundary curves:", fat_boundaries)
bone_muscle_interface = gmsh.model.addPhysicalGroup(1,
[tag[1] for tag in bone_boundaries],
name="Bone_Muscle")
muscle_fat = list(set(muscle_boundaries) - set(bone_boundaries))
muscle_fat_interface = gmsh.model.addPhysicalGroup(1,
[tag[1] for tag in muscle_fat],
name="Muscle_Fat")
fat_skin = list(set(fat_boundaries) - set(muscle_boundaries))
fat_skin_interface = gmsh.model.addPhysicalGroup(1,
[tag[1] for tag in fat_skin],
name="Fat_Skin")
skin_electrode = list(set(skin_boundaries) - set(fat_boundaries))
skin_electrode_interface = gmsh.model.addPhysicalGroup(1,
[tag[1] for tag in skin_electrode],
name="Skin_Electrode")
gmsh.model.occ.synchronize()
print("Muscle_Fat:", muscle_fat)
print("Fat_Skin:", fat_skin)
print("Skin_Electrode:", skin_electrode)
rect = gmsh.model.occ.addRectangle(-L0 / 2, -(r + t), 0, L0, t)
conductor_array = [gmsh.model.occ.addRectangle(IED / 2 + (Ed + IED) * i - (L0 / 2),
t - Et - (r + t), 0, Ed, Et) for i in range(count)]
gmsh.model.occ.synchronize()
print("Conductor array:", conductor_array)
# Subtracting to get embedded conductors
encaps_matrix = gmsh.model.occ.cut([(2, rect)],
[(2, tag) for tag in conductor_array],
removeObject=True, removeTool=False)
gmsh.model.occ.synchronize()
# Creating physical groups for surfaces
substrate_surfaces = [ent[1] for ent in encaps_matrix[0]]
substrate_tag = gmsh.model.addPhysicalGroup(2, substrate_surfaces,
name="Substrate")
conductor_tag = gmsh.model.addPhysicalGroup(2, conductor_array,
name="Conductor_Array")
gmsh.model.occ.synchronize()
boundaries = gmsh.model.getBoundary([(2, tag) for tag in substrate_surfaces],
oriented=False)
left_curves = []
right_curves = []
tol = 1e-6 # Tolerance
for dim, tag in boundaries:
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(dim,tag)
if abs(xmin + (L0 / 2)) < tol and abs(xmax + (L0 / 2)) < tol:
left_curves.append(tag)
elif abs(xmin - (L0 / 2)) < tol and abs(xmax - (L0 / 2)) < tol:
right_curves.append(tag)
if left_curves:
gmsh.model.addPhysicalGroup(1, left_curves, name="Left")
if right_curves:
gmsh.model.addPhysicalGroup(1, right_curves, name="Right")
print("Left curves:", left_curves)
print("Right curves:", right_curves)
gmsh.model.occ.synchronize()
# Group all band entities
#band_entities = encaps_matrix[0] + [(2, tag) for tag in conductor_array]
#gmsh.model.occ.translate(band_entities, -L0 / 2, -(r + t), 0)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(2)
gmsh.write("2d_forearm_electrode.msh")
gmsh.fltk.run()
gmsh.finalize()
And the next one attempts to move the band around the circle
import importlib.util
if importlib.util.find_spec("petsc4py") is not None:
import dolfinx
if not dolfinx.has_petsc:
print("This demo requires DOLFINx to be compiled with PETSc enabled.")
exit(0)
from petsc4py.PETSc import ScalarType # type: ignore
else:
print("This demo requires petsc4py")
exit(0)
from mpi4py import MPI
import numpy as np
import ufl
from dolfinx import fem, plot, default_scalar_type
from dolfinx.fem import form
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem.petsc import LinearProblem
import pyvista
msh_file = "2d_forearm_electrode.msh"
comm = MPI.COMM_WORLD
# (mm, Pa,)
L0 = 2 * np.pi * 30 # Band length
t = 1 # Band thickness
r = 45 # Radius of forearm
ν_eco = 0.49
E_eco = 68947.6 # Modulus
μ_eco = E_eco / (2.0 * (1.0 + ν_eco))
λ_eco = E_eco * ν_eco / ((1.0 + ν_eco) * (1.0 - 2.0 * ν_eco))
ν_cnt = 0.34
E_cnt = 1.1E+12
μ_cnt = E_cnt / (2.0 * (1.0 + ν_cnt))
λ_cnt = E_cnt * ν_cnt / ((1.0 + ν_cnt) * (1.0 - 2.0 * ν_cnt))
ν_rigid = 0.04
E_rigid = 1.1E+20
μ_rigid = E_rigid / (2.0 * (1.0 + ν_rigid))
λ_rigid = E_rigid * ν_rigid / ((1.0 + ν_rigid) * (1.0 - 2.0 * ν_rigid))
msh, cell_tags, facet_tags = gmshio.read_from_msh(msh_file, comm, 0, gdim=2)
tdim = msh.topology.dim
fdim = tdim - 1
# Create arrays matching number of cells
λ = np.zeros(msh.topology.index_map(tdim).size_local,
dtype=default_scalar_type)
μ = np.zeros_like(λ)
for cell_id, tag in enumerate(cell_tags.values):
if tag == 1:
λ[cell_id] = λ_eco
μ[cell_id] = μ_eco
elif tag == 2:
λ[cell_id] = λ_cnt
μ[cell_id] = μ_cnt
else: # All other forearm tissues rigid
λ[cell_id] = λ_rigid
μ[cell_id] = μ_rigid
V0 = fem.functionspace(msh, ("DG", 0))
λ_func = fem.Function(V0)
μ_func = fem.Function(V0)
λ_func.x.array[:] = λ
μ_func.x.array[:] = μ
print("Values: ", facet_tags.values)
print("Facet dim:", facet_tags.dim)
print(f"Whole face tag of type {facet_tags.name}=>", dir(facet_tags))
print(f"Geometry dimension is {msh.geometry.dim}")
V = fem.functionspace(msh, ("Lagrange", 2, (msh.geometry.dim, )))
left_facets = facet_tags.find(11)
right_facets = facet_tags.find(12)
left_dofs = fem.locate_dofs_topological(V, fdim, left_facets)
right_dofs = fem.locate_dofs_topological(V, fdim, right_facets)
u_left = np.array([L0 / 2, 2 * r], dtype=default_scalar_type)
u_right = np.array([-L0 / 2, 2 * r], dtype=default_scalar_type)
bc_left = fem.dirichletbc(u_left, left_dofs, V)
bc_right = fem.dirichletbc(u_right, right_dofs, V)
bcs = [bc_left, bc_right]
ds = ufl.Measure("ds", domain=msh)
def σ(u):
"""Return an expression for the stress σ given a displacement field"""
return (2.0 * μ_func * ufl.sym(ufl.grad(u))
+ λ_func * ufl.nabla_div(u) * ufl.Identity(len(u)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(σ(u), ufl.grad(v)) * ufl.dx
T = fem.Constant(msh, default_scalar_type((0, 0)))
f = fem.Constant(msh, default_scalar_type((0, 0)))
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
problem = LinearProblem(
a,
L,
bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
)
uh = problem.solve()
# Create the VTK mesh from your function space
topology, cell_types, geometry = plot.vtk_mesh(V)
# Create an unstructured grid (2D domain)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
# Attach displacement data (reshape correctly for 2D)
uh_array = uh.x.array.reshape((geometry.shape[0], msh.geometry.dim))
# If 2D, pad to 3 components (z = 0) because VTK expects 3-component vectors
if msh.geometry.dim == 2:
zeros = np.zeros((uh_array.shape[0], 1), dtype=uh_array.dtype)
uh_3 = np.hstack([uh_array, zeros])
else:
uh_3 = uh_array
# Ensure dtype is float64 for PyVista
uh_3 = uh_3.astype(np.float64)
grid["u"] = uh_3
"""uh_x = np.zeros(uh_3.shape)
uh_y = np.zeros(uh_3.shape)
uh_x[:, 0] = uh_3[:, 0]
uh_y[:, 1] = uh_3[:, 1]
grid["u_xy"] = uh_x + uh_y"""
grid.set_active_vectors("u")
# Warp geometry by displacement (for visualization)
# Factor scales how much deformation is shown (visual only)
warped = grid.warp_by_vector("u", factor=1.0)
# Initialize the plotter
p = pyvista.Plotter()
#p.add_mesh(grid, style="wireframe", color="black", line_width=0.5, label="Original")
##p.add_mesh(warped, scalars=np.linalg.norm(uh_array, axis=1),
## cmap="viridis", show_edges=True, label="Deformed")
#p.add_mesh(grid, scalars="u_xy", cmap="coolwarm", show_edges=True)
p.add_mesh(grid, style="wireframe", color="k")
p.add_mesh(warped, show_edges=True)
p.add_scalar_bar(title="u_xy (signed displacement)")
p.show_axes()
p.view_xy()
if not pyvista.OFF_SCREEN:
p.show()
else:
figure_as_array = p.screenshot("deformation.png")
In the end, there seem to be not movement on the band