Hello !
I’m relatively new to using DOLFINx and am currently working on a project involving the modeling of a lattice structure using Timoshenko beam elements. I’ve been trying to simulate the displacement of the structure under a compressive force applied to the top surface. However, I’m encountering an issue where the displacement results in infinity, and I’m unsure where I might be going wrong.Here’s a brief overview of what I’ve done so far:
I have successfully read the vertices and edges from an input file and created a mesh for the lattice structure.
I’ve defined the function space for displacement and rotation using appropriate finite element types.
The stiffness matrix (K) and load vector have been formulated using variational forms, taking into account both bending and shear deformations.
I’ve applied Dirichlet boundary conditions to fix the bottom nodes of the structure.
A compressive force of 100N has been applied to the top nodes, and I’ve attempted to solve the linear system to find the displacement field.
Despite these efforts, the displacement values are resulting in infinity. I suspect there might be an issue with the application of point forces or the assembly of the system, but I’m not sure how to resolve it. I would appreciate any advice on how to solve this issue .
import numpy as np
from mpi4py import MPI
import dolfinx
import basix
from dolfinx.mesh import create_mesh, locate_entities_boundary, create_interval
from dolfinx.fem import functionspace, Function, Constant, dirichletbc, locate_dofs_topological, form, apply_lifting, set_bc
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.io import XDMFFile
from petsc4py import PETSc
import ufl
from ufl import TrialFunction, TestFunction, dot, inner, grad, dx, split
# Define the vertices of the cube
vertices = np.array([
[-1, -1, -1],
[1, -1, -1],
[1, 1, -1],
[-1, 1, -1],
[-1, -1, 1],
[1, -1, 1],
[1, 1, 1],
[-1, 1, 1]
])
# Define the edges of the cube
edges = [
(0, 1),
(1, 2),
(2, 3),
(3, 0),
(4, 5),
(5, 6),
(6, 7),
(7, 4),
(0, 4),
(1, 5),
(2, 6),
(3, 7)
]
# Define the cell type and degree, including the shape
cell_type = basix.ufl.element("Lagrange", "interval", degree=1, shape=(3,))
# Create the mesh
domain = create_mesh(MPI.COMM_WORLD, edges, vertices, cell_type)
# Verify the mesh
print("Mesh created successfully.")
# Define the function space for displacement (w) and rotation (theta)
U = basix.ufl.element("Lagrange", domain.topology.cell_name(), 2) # Displacement
T = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1) # Rotation
V = functionspace(domain, basix.ufl.mixed_element([U, T]))
print("Function space defined successfully.")
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
(w, theta) = split(u)
(dw, dtheta) = split(v)
print("Trial and test functions defined successfully.")
# Material properties
E = 70e3 # Young's modulus in MPa
nu = 0.3 # Poisson's ratio
G = E / (2 * (1 + nu)) # Shear modulus in MPa
r = 0.05 # Radius of the circular cross-section
A = np.pi * r**2 # Cross-sectional area
I = np.pi * r**4 / 4 # Second moment of area
kappa = 5 / 6 # Shear correction factor
# Define the stiffness matrix (K) bilinear form
k_form = (E * I * inner(grad(theta), grad(dtheta)) * dx +
kappa * G * dot(grad(w)[0] - theta, grad(dw)[0] - dtheta) * dx)
print("Stiffness matrix bilinear form defined successfully.")
# Define linear form for assembling the stiffness matrix
constant = Constant(domain, PETSc.ScalarType(1.0))
l_form = constant * w * dx
print(" linear form defined successfully.")
# Locate boundary vertices based on minimum z-coordinate
z_min = np.min(vertices[:, 2])
boundary_facets = locate_entities_boundary(domain, domain.topology.dim - 1, lambda x: np.isclose(x[2], z_min))
# Apply Dirichlet boundary condition (e.g., w = 0 on the boundary)
bc = dirichletbc(PETSc.ScalarType(0), locate_dofs_topological(V.sub(0), domain.topology.dim - 1, boundary_facets), V.sub(0))
# Assemble the stiffness matrix (K)
a = form(k_form)
K = assemble_matrix(a, bcs=[bc])
K.assemble()
print("Stiffness matrix assembled successfully.")
# Locate the top nodes based on maximum z-coordinate
z_max = np.max(vertices[:, 2])
top_nodes = locate_entities_boundary(domain, domain.topology.dim - 1, lambda x: np.isclose(x[2], z_max))
# Map the top nodes to DOFs
top_dofs = locate_dofs_topological(V.sub(0), domain.topology.dim - 1, top_nodes)
# Define the force vector (f) with a point load of 100N applied to the top nodes
f = Function(V)
f.vector.set(0.0)
for dof in top_dofs:
f.vector.array[dof] += 100.0
# Apply boundary conditions to the force vector
apply_lifting(f.vector, [a], bcs=[[bc]], scale=1.0)
f.vector.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(f.vector, [bc])
# Solve the linear system K * u = f
u = Function(V)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(K) # K is already a PETSc Mat
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.solve(f.vector, u.vector)
u.vector.ghostUpdate()
print("Displacement solution obtained successfully.")
# Print the results
print("Displacement solution:")
for i, displacement in enumerate(u.vector.array):
print(f"Node {i}: {displacement}")
