Hi, is there a way of applying a force load for 3D linear elasticity if the node number (global) is known?
I am trying to modify the minimal 3D linear elasticity code below. I took out the self weight term in the variational formulation.
I remember there was a ‘pointwise’ option for Diritchlet BCs in the legacy FEnics, but it also required the use of a marker function to find the node. Is there a function I can use to apply the force at a node if the global node number is known? Maybe accessing the force vector of the linear system?
# Common libraries
import numpy as np
# Fenicsx FEA libraries
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
mu = 1
lambda_ = 1.25
domain = mesh.create_unit_cube(MPI.COMM_WORLD, nx=1, ny=1, nz=1, cell_type=mesh.CellType.tetrahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
def clamped_boundary(x):
return np.isclose(x[0], 0)
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)
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(T, v) * ds
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
uh.name = "Deformation"
xdmf.write_function(uh)
