I want to see the force-displacement curve for the points at centre of my mesh.
I thought about doing something like this, but I don’t know how to go from there:
def middle_region(x): # region at middle-center
halflengthmiddlearea = length/50
return (x[0] >= length/2 - halflengthmiddlearea) & (x[0] <= length/2 + halflengthmiddlearea) & (
x[1] >= height/2 - halflengthmiddlearea) & (x[1] <= height/2 + halflengthmiddlearea)
def left_region(x): # region at left-center
halflengthmiddlearea = length/50
return (x[0] >= 0) & (x[0] <= halflengthmiddlearea) & (
x[1] >= height/2 - halflengthmiddlearea) & (x[1] <= height/2 + halflengthmiddlearea)
cells_centre = dolfinx.mesh.locate_entities(domain, domain.topology.dim, middle_region)
cells_left = dolfinx.mesh.locate_entities(domain, domain.topology.dim, left_region)
Also, i tried using mesh.locate_entities (inspired from Defining subdomains for different materials — FEniCSx tutorial) but i think i am doing something wrong
# evaluate displacement at points of interest
for point in points_of_interest:
cell_index = mesh.locate_entities(domain, dim, point)
displacement = u_h.eval(point, cell_index)
displacements.append(displacement)
this is my full code:
import numpy as np # type: ignore
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction # type: ignore
from mpi4py import MPI # type: ignore
from dolfinx import fem, io, mesh
from dolfinx.mesh import create_rectangle, CellType
import dolfinx.fem.petsc
import ufl # type: ignore
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from petsc4py import PETSc # type: ignore
from mpi4py import MPI # type: ignore
# Define temporal parameters
t = 0 # Start time
T = 5.0 # Final time
timesteps = 100
dt = T / timesteps # time step size
# Domain parameters
length, height = 10, 1.0
Nx, Ny = 50, 10
# Create mesh
domain = create_rectangle(
MPI.COMM_WORLD, # parallel communicator
[np.array([0, 0]), np.array([length, height])],
[Nx, Ny],
cell_type=CellType.quadrilateral,
)
dim = domain.topology.dim
print(f"Mesh topology dimension d={dim}.")
degree = 2
shape = (dim,) # this means we want a vector field of size `dim`
# Define the function space
V = fem.functionspace(domain, ("P", degree, shape))
u = TrialFunction(V)
v = TestFunction(V)
u_h = fem.Function(V, name="Displacement")
u_n = fem.Function(V)
def middle_region(x):
halflengthmiddlearea = length/50
return (x[0] >= length/2 - halflengthmiddlearea) & (x[0] <= length/2 + halflengthmiddlearea) & (
x[1] >= height/2 - halflengthmiddlearea) & (x[1] <= height/2 + halflengthmiddlearea)
def left_region(x):
halflengthmiddlearea = length/10
return (x[0] >= 0) & (x[0] <= halflengthmiddlearea) & (
x[1] >= height/2 - halflengthmiddlearea) & (x[1] <= height/2 + halflengthmiddlearea)
cells_centre = dolfinx.mesh.locate_entities(domain, domain.topology.dim, middle_region)
cells_left = dolfinx.mesh.locate_entities(domain, domain.topology.dim, left_region)
# Lamé parameters (for a given material)
E = 10.0 # Young's modulus
nu = 0.3 # Poisson's ratio
lambda_ = (E * nu) / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
# Define strain and stress tensors
def epsilon(u):
return sym(grad(u))
def sigma(u):
return lambda_ * tr(epsilon(u)) * Identity(len(u)) + 2 * mu * epsilon(u)
# Dirichlet boundary conditions or None
bcs = []
# Mark the right boundary with ID 1, here the force will pull
fdim = domain.topology.dim - 1
right_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.isclose(x[0], length))
right_marker = mesh.meshtags(domain, fdim, right_facets, 1)
# Define the boundary measure using the right marker
ds = Measure("ds", domain=domain, subdomain_data=right_marker)
# Locate top and bottom facets for friction
top_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_and(np.isclose(x[1], height), x[0] <= length / 2))
bottom_facets = mesh.locate_entities_boundary(domain, fdim, lambda x: np.logical_and(np.isclose(x[1], 0), x[0] <= length / 2))
# Combine top and bottom markers into a single measure
friction_marker = mesh.meshtags(domain, fdim, np.concatenate([top_facets, bottom_facets]), 1)
ds_friction = Measure("ds", domain=domain, subdomain_data=friction_marker)
# Weak formulation
f_x_value = 1.0 # magnitude of the force in the x-direction
friction_coefficient = 0.5
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = fem.Constant(domain, np.array([f_x_value, 0.0]))
dx = Measure("dx", domain=domain)
# linear form
a = (
inner(sigma(u), epsilon(v)) * dx +
(1 / dt) * inner(u, v) * dx +
friction_coefficient * inner(u, v) * ds_friction(1)
)
# billinear form
L = inner(f, v) * ds(1) + (1 / dt) * inner(u_n, v) * dx
bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = assemble_matrix(bilinear_form, bcs=bcs)
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
# define the problem
problem = fem.petsc.LinearProblem(
a, L, u=u_h, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
vtk = io.VTKFile(domain.comm, "results/dynamic.pvd", "w")
displacements = []
# points along the horizontal centre line of the domain
points_of_interest = [[x, height / 2, 0] for x in np.linspace(0, length, 10)]
for i in range(timesteps):
print(f"Time step {i}/{timesteps}")
t += dt
if i == 30:
f.value[0] *= 2
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
# Solve linear problem
solver.solve(b, u_h.x.petsc_vec)
u_h.x.scatter_forward()
# Update solution at previous time step (u_n)
u_n.x.array[:] = u_h.x.array
# evaluate displacement at points of interest
for point in points_of_interest:
cell_index = mesh.locate_entities(domain, dim, middle_region)
displacement = u_h.eval(point, cell_index)
displacements.append(displacement)
# Write solution to file
vtk.write_function(u_h, t)
vtk.close()
I hope I made it understandable what I am struggling with and that someone might be able to help 