‘’’
dolfinx: 0.8.0
‘’’
import numpy as np
# Define the 3x3 force and rotation matrices
K_fx = np.array([[K_fxfx, K_fxfy, K_fxfz],
[K_fyfx, K_fyfy, K_fyfz],
[K_fzfx, K_fzfy, K_fzfz]])
K_mx = np.array([[K_mxfx, K_mxfy, K_mxfz],
[K_myfx, K_myfy, K_myfz],
[K_mzfx, K_mzfy, K_mzfz]])
# Initialize a 6x6 zero matrix
K_combined = np.zeros((6, 6))
# Place the 3x3 matrices in the appropriate positions
K_combined[:3, :3] = K_fx
K_combined[3:, 3:] = K_mx
print(K_combined)
\[
\begin{bmatrix}
K_{fxfx} & K_{fxfy} & K_{fxfz} & K_{fxmx} & K_{fxmy} & K_{fxmz} \\
K_{fyfx} & K_{fyfy} & K_{fyfz} & K_{fymx} & K_{fymy} & K_{fymz} \\
K_{fzfx} & K_{fzfy} & K_{fzfz} & K_{fzmx} & K_{fzmy} & K_{fzmz} \\
K_{mxfx} & K_{mxfy} & K_{mxfz} & K_{mxmx} & K_{mxmy} & K_{mxmz} \\
K_{myfx} & K_{myfy} & K_{myfz} & K_{mymx} & K_{mymy} & K_{mymz} \\
K_{mzfx} & K_{mzfy} & K_{mzfz} & K_{mzmx} & K_{mzmy} & K_{mzmz}
\end{bmatrix}
\]
So then, I have this formula, which is supposed to be a description of how to combine the forces and moments compartments, and I am really not sure if this is the correct way yet of doing things however it does look like it has a chance to work so far in my assesment…
Is it possible to do this type of formulation with dolfinx?
So far I made an attempt at doing a similar type of thing but there seems to be a difficulty at solution time. I’m not sure what went wrong here yet, so far the ufl system is still a bit mysterious to me and I have the gist of what those computations mean but there must be something thats causing a congruency and or shape issue at solution time?
Could I get some help with this?
import dolfinx
from dolfinx import mesh, fem
from mpi4py import MPI
import basix
import numpy as np
from petsc4py import PETSc
import ufl
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import dirichletbc, locate_dofs_geometrical
L = 72
W = 4.0
# Define the mesh with tetrahedral cells
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
[20, 6, 6], cell_type=mesh.CellType.tetrahedron)
# Define elements for displacements and rotations
element_disp = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
element_rot = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(3,))
# Create function spaces for displacements and rotations
V_disp = fem.functionspace(domain, element_disp)
V_rot = fem.functionspace(domain, element_rot)
# Create functions to hold the force vectors at each node
force_disp = fem.Function(V_disp)
force_rot = fem.Function(V_rot)
# Initialize the force vectors to zero
force_disp.vector.set(0.3)
force_rot.vector.set(0.0)
# Define functions for displacement and rotation
u_disp = fem.Function(V_disp)
u_rot = fem.Function(V_rot)
v_disp = fem.Function(V_disp)
v_rot = fem.Function(V_rot)
# Define UFL expressions for the displacement and rotation components
grad_u_disp = ufl.grad(u_disp) # Gradient of displacement u
grad_v_disp = ufl.grad(v_disp) # Gradient of displacement v
# Symmetric part of the gradient (strain tensor) for displacements
epsilon_u = ufl.sym(grad_u_disp)
epsilon_v = ufl.sym(grad_v_disp)
# Bilinear form for displacements
a_disp = fem.form(ufl.inner(epsilon_u, epsilon_v) * ufl.dx)
# (Example) Bilinear form for rotations, assuming small rotations
a_rot = fem.form(ufl.inner(u_rot, v_rot) * ufl.dx)
# Total bilinear form (sum of displacement and rotation contributions)
a = fem.form(ufl.inner(epsilon_u, epsilon_v) * ufl.dx + ufl.inner(u_rot, v_rot) * ufl.dx)
# Define the linear form (force vector)
L_disp = fem.form(ufl.dot(force_disp, v_disp) * ufl.dx)
L_rot = fem.form(ufl.dot(force_rot, v_rot) * ufl.dx)
# Define the total linear form
L = fem.form(ufl.dot(force_disp, v_disp) * ufl.dx + ufl.dot(force_rot, v_rot) * ufl.dx)
# Define the boundary condition
def boundary(x):
return np.isclose(x[0], 0)
# Define the Dirichlet boundary value for vector-valued spaces
zero_disp = fem.Constant(domain, PETSc.ScalarType((0.0, 0.0, 0.0)))
zero_rot = fem.Constant(domain, PETSc.ScalarType((0.0, 0.0, 0.0)))
# Apply Dirichlet boundary conditions to the displacements and rotations
bc_disp = dirichletbc(zero_disp, locate_dofs_geometrical(V_disp, boundary), V_disp)
bc_rot = dirichletbc(zero_rot, locate_dofs_geometrical(V_rot, boundary), V_rot)
# Combine boundary conditions
bcs = [bc_disp, bc_rot]
# Define the linear problem
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# Solve the linear problem
u_sol = problem.solve()
# Split the solution into displacement and rotation components
u_disp, u_rot = u_sol.split()
# Output the results
print(f"Displacement:\n{u_disp.x.array}")
print(f"Rotation:\n{u_rot.x.array}")
print(f'dolfinx version: {dolfinx.__version__}')