How to apply point forces by modifying the assembled matrix for static linear elasticity

Hi Fenicsx community,

I’m trying to modify the assembled matrix from the linear part in order to apply a point force. I saw from here : Dirac delta distribution (dolfinx) that the Dirac delta distribution is not yet implemented in fenicsx, unless I misunderstood.

But like said in https://fenicsproject.org/qa/14258/how-to-apply-pointsource-to-linearvariationalproblem/, there is another option, add the force on assembled vectors.

So far, I’ve tried to do the second solution, it works, but only if I solve the problem with numpy, by transferring the dense matrix to python, but it is not an efficient method. Here is what I did :

``````from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import dx, ds, Measure, TrialFunction, TestFunction, sym, grad, div, Identity, inner, dot
from dolfinx import mesh, fem, plot, io
import numpy as np

domain = mesh.create_rectangle(MPI.COMM_SELF, ((-2, 0), (2, 1)), (50,30), mesh.CellType.quadrilateral)
V2 = fem.VectorFunctionSpace(domain, ("CG", 2))

E = fem.Constant(domain, ScalarType(1))
nu = fem.Constant(domain, ScalarType(0.3))
lamda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / (2 * (1 + nu))
T = fem.Constant(domain, ScalarType((0, 0)))

boundary_dofs = fem.locate_dofs_topological(V2, f_dim, lambda x: np.isclose(x[0], -2))
bc = fem.dirichletbc(np.zeros([2], dtype=ScalarType), boundary_dofs, V2)

dofs = fem.locate_dofs_geometrical(V2, lambda x: np.isclose(x[0], 2) & np.isclose(x[1], 0.5))

u = TrialFunction(V2)
v = TestFunction(V2)

def eps(u):

def sigma(u):
return (lamda * div(u) * Identity(2) + 2 * mu * eps(u))

a = inner(sigma(u), eps(v)) * dx
L = dot(T, v) * ds

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.b.setValue(dofs, -1)  # Assign the force on the good dof

uh = problem.solve()
print((uh.x.array == 0).all())  # return True

print(problem.b.getArray()[dofs])  # After solving the problem it is = 0 instead of -1 :(
``````

`solve` calls assemble vector which would overwrite any entries in vector:

You would need to use the PETSc KSP interface directly, or make a modified version of the solve command, that applies this entry to the vector after assembling and applying lifting.

1 Like

Thank you very much @dokken it works perfectly.

Here the code that I use :

``````class MyLinearProblem(fem.petsc.LinearProblem):
def __init__(self, *args, **kwargs):
if 'point_force_dofs' in kwargs and 'point_force_values' in kwargs:
point_force_dofs = kwargs['point_force_dofs']
point_force_values = kwargs['point_force_values']
if not isinstance(type(point_force_dofs), np.ndarray):
point_force_dofs = np.array([point_force_dofs], dtype='int64')
if not isinstance(type(point_force_values), np.ndarray):
point_force_values = np.array([point_force_values])
point_force_dofs = point_force_dofs.flatten()
point_force_values = point_force_values.flatten()

if len(point_force_dofs) != len(point_force_values):
raise f"The length of point_force_dofs {len(point_force_values)} is not equal to the length of " \
f"point_force_values"

del kwargs['point_force_dofs'], kwargs['point_force_values']
else:
point_force_dofs = None
point_force_values = None
self._point_force_dofs = point_force_dofs
self._point_force_values = point_force_values

super().__init__(*args, **kwargs)

def solve(self):
# Assemble lhs
self._A.zeroEntries()
from dolfinx.fem.petsc import _assemble_matrix_mat
_assemble_matrix_mat(self._A, self._a, bcs=self.bcs)
self._A.assemble()

# Assemble rhs
with self._b.localForm() as b_loc:
b_loc.set(0)
from dolfinx.fem.petsc import assemble_vector
assemble_vector(self._b, self._L)

# Apply boundary conditions to the rhs
from dolfinx.fem.petsc import apply_lifting
apply_lifting(self._b, [self._a], bcs=[self.bcs])
from dolfinx.fem.petsc import set_bc
set_bc(self._b, self.bcs)

if self._point_force_values is not None and self._point_force_dofs is not None:
self._b.setValue(self._point_force_dofs, self._point_force_values)

# Solve linear system and update ghost values in the solution
self._solver.solve(self._b, self._x)
self.u.x.scatter_forward()

return self.u
``````

Have a nice evening !

1 Like