dolfinx 0.8.0…
Somewhere on this forum I read that there should be a way with the latest build of dolfinx to apply a point load via modification of the Right Hand Side Vector, and that the methodology to do so has changed since fenicsx…
This is my first attempt to set up such a scenario:
import numpy as np
import pandas as pd
import pyvista
from mpi4py import MPI
from petsc4py import PETSc
import ufl
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem import petsc
import dolfinx
print(f'dolfinx version:{dolfinx.__version__}')
# --------------------------------------------------------------------
# Problem parameters and mesh
# --------------------------------------------------------------------
# Define dimensions in inches
L_domain = 120.0 # length of the beam
W = 4.0 # width
H = 8.0 # height
# Number of elements (here chosen so that elements are ~0.8" on a side)
nx, ny, nz = 150 * 2, 5 * 2, 10 * 2
degree = 1 # polynomial degree for the function space
# Create a box mesh (a hexahedral mesh)
domain = mesh.create_box(MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([L_domain, W, H])],
[nx, ny, nz],
cell_type=mesh.CellType.hexahedron)
# Define a vector-valued Lagrange finite element space
V = fem.functionspace(domain, ("Lagrange", degree, (domain.geometry.dim, )))
# --------------------------------------------------------------------
# Boundary conditions (pinned at both ends)
# --------------------------------------------------------------------
def clamped_boundary(x):
return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], L_domain))
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D,
fem.locate_dofs_topological(V, fdim, boundary_facets),
V)
# --------------------------------------------------------------------
# Material properties and variational formulation
# --------------------------------------------------------------------
# Material properties for structural wood (psi)
E = 2e6 # Young's modulus
nu = 0.3 # Poisson's ratio
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
# Define strain and stress functions
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Density of structural wood (lb/in^3)
rho = 0.01736
# Gravitational body force (lb/in^3)
f_z_gravity = -rho
# Distributed load: given as lb/inch along the beam length.
w = -0.5 # lb/inch
total_load = w * L_domain # total load (lb)
volume = L_domain * W * H
f_z_load = total_load / volume
# Total body force in z-direction (lb/in^3)
f_z_total = f_z_load + f_z_gravity
# Define the body force vector as a constant
f = fem.Constant(domain, default_scalar_type((0, 0, f_z_total)))
# Traction is zero on the Neumann part of the boundary (if any)
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
# Define the variational forms using fem.form
a = fem.form(ufl.inner(sigma(u), epsilon(v)) * ufl.dx)
L = fem.form(ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds)
# Create the stiffness matrix
A = fem.petsc.create_matrix(a)
# Assemble the stiffness matrix with boundary conditions
fem.petsc.assemble_matrix(A, a, bcs=[bc])
A.assemble()
# Create the right-hand side vector
b = fem.petsc.create_vector(L)
# Assemble the right-hand side vector
fem.petsc.assemble_vector(b, L)
# Apply lifting for the boundary conditions
fem.apply_lifting(b, [a], [bc])
# Update ghost values
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Set the boundary conditions on the vector
fem.set_bc(b, [bc])
# --------------------------------------------------------------------
# Modify RHS to include an additional point load
# --------------------------------------------------------------------
# Choose the point at which to apply the load (here, at the midpoint)
point = np.array([L_domain / 2, W / 2, H / 2], dtype=default_scalar_type)
# Get the coordinates of the degrees of freedom (reshape since V is vector-valued)
dof_coords = V.tabulate_dof_coordinates().reshape((-1, 3))
# Find the index of the node closest to the desired point
node_index = np.argmin(np.linalg.norm(dof_coords - point, axis=1))
# For a 3D vector-valued space the dofs are typically arranged per node.
# Here, we assume that the ordering is [u_x, u_y, u_z] at each node.
# Hence, the global index for the z-component at the chosen node is:
load_dof = 3 * node_index + 2
# Define the point load magnitude (e.g., -10 lb in the z-direction)
P = -10.0
# Add the point load to the assembled RHS vector. In a distributed setting,
# setValue will only affect the process that owns that dof.
b.setValue(load_dof, P, addv=PETSc.InsertMode.ADD)
b.assemble()
# --------------------------------------------------------------------
# Solve the linear system A u = b
# --------------------------------------------------------------------
uh = fem.Function(V)
# Create a PETSc KSP solver and configure it for a direct LU solve
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType("preonly")
pc = solver.getPC()
pc.setType("lu")
# Solve the system
solver.solve(b, uh.vector)
uh.x.scatter_forward() # update ghost values if needed
print("Minimum displacement:", uh.x.array.min())
# --------------------------------------------------------------------
# Post-processing: Export results and visualization
# --------------------------------------------------------------------
# ~~~~
dolfinx version:0.8.0
Traceback (most recent call last):
File "/home/prusso/dolfinx-testcase1-beamsimple_w/dolfinx_testcase2.py", line 106, in <module>
fem.apply_lifting(b, [a], [bc])
File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/assemble.py", line 347, in apply_lifting
_bcs = [[bc._cpp_object for bc in bcs0] for bcs0 in bcs]
File "/home/prusso/spack/var/spack/environments/fenicsx-env/.spack-env/view/lib/python3.10/site-packages/dolfinx/fem/assemble.py", line 347, in <listcomp>
_bcs = [[bc._cpp_object for bc in bcs0] for bcs0 in bcs]
TypeError: 'DirichletBC' object is not iterable
So I have this difficulty to work through as shown in the above traceback. Also as far as the point load I want to apply the point load in the downward direction. So far the comment implies the Z-direction which I assume so far is along the length of the beam and the -Y-direction should be downward as far as I know. So far I would like to apply the point load in addition to the distributed load f_z_load which is working quite well in a more isolated scripted environment…
Could I get some help to get this point load applied with displacements coming out to look at?