Hi, I am trying to apply fem analysis on a 1d nonlinear timoshenko spatial beam, where there is 9 dofs. I tried to implement in the following way
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import pyvista
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, io, mesh
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner, grad)
from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 1.0
points = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]])
# Create a 1D mesh for the beam
beam_mesh = mesh.create_interval(MPI.COMM_WORLD, 5, points[:,0])
displacement1 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
displacement2 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
displacement3 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
curvature1 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
curvature2 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
curvature3 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
rotation1 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
rotation2 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
rotation3 = element("Lagrange", beam_mesh.basix_cell(), degree=1)
element = mixed_element([displacement1, displacement2, displacement3, curvature1, curvature2, curvature3, rotation1, rotation2, rotation3])
#defining functionspace
V = fem.functionspace(beam_mesh, element)
#define functions for displacements
#u = fem.Function(V)
(z1, z2, z3, k1, k2, k3, p1, p2, p3) = TrialFunctions(V)
(u1, u2, u3, v1, v2, v3, w1, w2, w3) = TestFunctions(V)
x = SpatialCoordinate(beam_mesh)
G = fem.Constant(beam_mesh, 2.50e9)
A = fem.Constant(beam_mesh, 1.0e-6)
K_s = fem.Constant(beam_mesh, 0.851)
E = fem.Constant(beam_mesh, 200e9)
P = fem.Constant(beam_mesh, -10.0)
G = fem.Constant(beam_mesh, 75e9)
dx = Measure("dx", beam_mesh)
#domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 0.001, 0.001]], [5, 1, 1], mesh.CellType.hexahedron)
#V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))
# several derivatives with respect to x, y, z
dz1_dx = grad(z1)[0]
du1_dx = grad(u1)[0]
strain_11 = dz1_dx + 0.5 * (dz1_dx)**2
strain_12 = 0.5 * p2
strain_13 = 0.5 * p3
del_strain11 = du1_dx + 0.5 * (du1_dx)**2
del_strain12 = 0.5 * w2
del_strain13 = 0.5 * w3
stress_11 = E * strain_11
stress_12 = G * strain_12
stress_13 = G * strain_13
a = inner(del_strain11, stress_11) * dx + inner(del_strain12, stress_12) * dx + inner(del_strain13, stress_13) * dx
#Next, we define the body force on the reference configuration (B), and nominal traction (T).
B = fem.Constant(beam_mesh, default_scalar_type((0, 0, 0)))
T = fem.Constant(beam_mesh, default_scalar_type((0, 0, 0)))
# Function returning True for points on the left boundary
def left_end(x):
return np.isclose(x[0], 0.0)
# 2. Locate the degrees of freedom (DoFs) at the left boundary for both displacement and rotation
fdim = beam_mesh.topology.dim - 1 # Geometrical dimension for the boundary (1D case: 0-dimensional points)
facets_left = mesh.locate_entities_boundary(beam_mesh, fdim, left_end)
# Locate DoFs for both displacement (w) and rotation (phi)
V_z1, _ = V.sub(0).collapse() # Subspace for displacement
V_z2, _ = V.sub(1).collapse() # Subspace for rotation
V_z3, _ = V.sub(2).collapse() # Subspace for displacement
V_k1, _ = V.sub(3).collapse() # Subspace for rotation
V_k2, _ = V.sub(3).collapse() # Subspace for displacement
V_k3, _ = V.sub(5).collapse() # Subspace for rotation
V_p1, _ = V.sub(6).collapse() # Subspace for displacement
V_p2, _ = V.sub(7).collapse() # Subspace for rotation
V_p3, _ = V.sub(8).collapse() # Subspace for displacement
dofs_z1_left = fem.locate_dofs_topological((V.sub(0), V_z1), fdim, facets_left)
dofs_z2_left = fem.locate_dofs_topological((V.sub(1), V_z2), fdim, facets_left)
dofs_z3_left = fem.locate_dofs_topological((V.sub(2), V_z3), fdim, facets_left)
dofs_k1_left = fem.locate_dofs_topological((V.sub(3), V_k1), fdim, facets_left)
dofs_k2_left = fem.locate_dofs_topological((V.sub(4), V_k2), fdim, facets_left)
dofs_k3_left = fem.locate_dofs_topological((V.sub(5), V_k3), fdim, facets_left)
dofs_p1_left = fem.locate_dofs_topological((V.sub(6), V_p1), fdim, facets_left)
dofs_p2_left = fem.locate_dofs_topological((V.sub(7), V_p2), fdim, facets_left)
dofs_p3_left = fem.locate_dofs_topological((V.sub(8), V_p3), fdim, facets_left)
# 3. Apply Dirichlet boundary conditions for z, k, p = 0
# Define zero-valued functions for the boundary conditions
zero_z1 = fem.Function(V_z1) # Function for displacement boundary condition (w = 0)
zero_z2 = fem.Function(V_z2) # Function for rotation boundary condition (phi = 0)
zero_z3 = fem.Function(V_z3)
zero_k1 = fem.Function(V_k1)
zero_k2 = fem.Function(V_k2)
zero_k3 = fem.Function(V_k3)
zero_p1 = fem.Function(V_p1)
zero_p2 = fem.Function(V_p2)
zero_p3 = fem.Function(V_p3)
zero_z1.interpolate(lambda x: np.zeros_like(x[0]))
zero_z2.interpolate(lambda x: np.zeros_like(x[0]))
zero_z3.interpolate(lambda x: np.zeros_like(x[0]))
zero_k1.interpolate(lambda x: np.zeros_like(x[0]))
zero_k2.interpolate(lambda x: np.zeros_like(x[0]))
zero_k3.interpolate(lambda x: np.zeros_like(x[0]))
zero_p1.interpolate(lambda x: np.zeros_like(x[0]))
zero_p2.interpolate(lambda x: np.zeros_like(x[0]))
zero_p3.interpolate(lambda x: np.zeros_like(x[0]))
# Create the boundary condition objects
bc_z1_left = fem.dirichletbc(zero_z1, dofs_z1_left, V.sub(0)) # BC for displacement
bc_z2_left = fem.dirichletbc(zero_z2, dofs_z2_left, V.sub(1)) # BC for rotation
bc_z3_left = fem.dirichletbc(zero_z3, dofs_z3_left, V.sub(2)) # BC for rotation
bc_k1_left = fem.dirichletbc(zero_k1, dofs_k1_left, V.sub(3)) # BC for displacement
bc_k2_left = fem.dirichletbc(zero_k2, dofs_k2_left, V.sub(4)) # BC for rotation
bc_k3_left = fem.dirichletbc(zero_k3, dofs_k3_left, V.sub(5)) # BC for rotation
bc_p1_left = fem.dirichletbc(zero_p1, dofs_p1_left, V.sub(6)) # BC for displacement
bc_p2_left = fem.dirichletbc(zero_p2, dofs_p2_left, V.sub(7)) # BC for rotation
bc_p3_left = fem.dirichletbc(zero_p3, dofs_p3_left, V.sub(8)) # BC for rotation
# Combine the boundary conditions
bcs = [bc_z1_left, bc_z2_left,bc_z3_left, bc_k1_left, bc_k2_left, bc_k3_left, bc_p1_left, bc_p2_left, bc_p3_left]
def right_end(x):
return np.isclose(x[0], L)
facets_right = mesh.locate_entities_boundary(beam_mesh, fdim, right_end)
# Check if any facets were found at the right end
if len(facets_right) == 0:
raise ValueError("No boundary facets found at the right end of the beam.")
# Create an array of markers (use 1 for the right end)
values_right = np.full(len(facets_right), 1, dtype=np.int32)
# Create the meshtags object for the boundary
boundary_markers = mesh.meshtags(beam_mesh, fdim, facets_right, values_right)
# Define boundary measure with subdomain markers
ds = Measure("ds", domain=beam_mesh, subdomain_data=boundary_markers)
# Define the body force term in the weak form
body_force_term = inner(u1, B[0]) * dx + inner(u2, B[1]) * dx + inner(u3, B[2]) * dx
# Define the nodal force term in the weak form (traction)
traction_term = inner(u1, T[0]) * ds(1) + inner(u2, T[1]) * ds(1) + inner(u3, T[2]) * ds(1)
# Complete the right-hand side of the weak form
rhs = body_force_term + traction_term
# Full variational form combining strain and body/nodal forces
F = a - rhs
u = fem.Function(V)
from dolfinx.jit import ffcx_jit
# Increase the timeout for JIT compilation
jit_options = {"timeout": 600} # Increase the timeout (in seconds) to 10 minutes
# Use the updated JIT options when creating the form
problem = NonlinearProblem(F, u, bcs, jit_options=jit_options)
after some times, it is showing errors like jit timeout.
what might be the issue here?