How to apply a force load at a known node number

[alex24]

Hi Jorgen, thank you for the prompt reply.

I’m trying to assemble the linear system so I can use the following:

How to apply a pointwise force using dolfinx

I figured out what the mistake was. The degrees-of-freedom of the point I want to find through “rd = dolfinx.fem.locate_dofs_geometrical(V=V, marker=lambda x: np.isclose(x[0], 2.0) & np.isclose(x[1], 0.0))” is actually The global index of the node. So if I want to apply a x-direction force on the node, I should use “b.setValue(2*rd, 1)” instead of “b.setValue(rd, 1)”. Still thanks for your reply.

Based on your code here:
https://jsdokken.com/dolfinx-tutorial/chapter2/heat_code.html

And the answer given in:

How to apply a pointwise force using dolfinx

Hi, everyone! I want to apply a force on a specify node of the mesh in only one direction (e.g x-direction) in linear elasticity problem. How should I do this?

I wrote the following code (below). However I get this solver error:
File ~/.local/lib/python3.10/site-packages/spyder_kernels/customize/utils.py:209 in exec_encapsulate_locals
exec_fun(compile(code_ast, filename, “exec”), globals)

File ~/files/fenics_point_loading_example.py:88
uh = solver.solve(b, u.x.petsc_vec)

File PETSc/KSP.pyx:395 in petsc4py.PETSc.KSP.solve

Error: error code 73
[0] KSPSolve() line 1085 in ./src/ksp/ksp/interface/itfunc.c
[0] KSPSolve_Private() line 850 in ./src/ksp/ksp/interface/itfunc.c
[0] KSPSetUp() line 406 in ./src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 1017 in ./src/ksp/pc/interface/precon.c
[0] PCSetUp_LU() line 87 in ./src/ksp/pc/impls/factor/lu/lu.c
[0] MatGetOrdering() line 177 in ./src/mat/order/sorder.c
[0] Object is in wrong state
[0] Not for unassembled matrix

from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, default_scalar_type, io
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
import numpy as np



mu = 1
lambda_ = 1.25

rho = 1
g = 10


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

f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds


a_compiled = fem.form(a)
L_compiled = fem.form(L)


# Assemble system, applying boundary conditions
A = assemble_matrix(a_compiled, bcs=[bc])
#A.assemble()

b = assemble_vector(L_compiled)
apply_lifting(b, [a_compiled], bcs=[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])

global_node = 5       # Example
b.setValue(3*global_node,   1)
b.setValue(3*global_node+1, 1)
b.setValue(3*global_node+2, 1)


# Create solution function
u = fem.Function(V)


# Solver
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

# Compute solution
uh = solver.solve(b, u.x.petsc_vec)



with io.XDMFFile(domain.comm, "deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)