Hi Everyone,

I think that I made a progress. By combining information from Dirac delta distribution (dolfinx) with Cahn-Hilliard equation demo I realized that u function in Dokkenās example is not TrialFunction as it is in many demos, but a separate function.

I attached a demo code below. Letās assume that you have a square 10x10 domain, and want to put two sources of intensities 1 and -1 on coordinates (2.5, 5) and (7.5, 5).

```
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (10.0, 10.0)), n=(100, 100),
cell_type=mesh.CellType.triangle,)
V = fem.FunctionSpace(msh, ("Lagrange", 1))
facets = mesh.locate_entities_boundary(msh, dim=1, marker=lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 10.0)),
np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 10.0))))
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = fem.Function(V)
dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x.T, [2.5, 5, 0]).all(axis=1))
f.x.array[dofs] = 1
dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x.T, [7.5, 5, 0]).all(axis=1))
f.x.array[dofs] = -1
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
try:
import pyvista
pyvista.start_xvfb()
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
grid.point_data["u"] = uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=False)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show(cpos='xy')
except ModuleNotFoundError:
print("'pyvista' is required to visualise the solution")
print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")
```

Result of modelling:

Please note that this method works only if the location of point sources aligns with the mesh.

Right now I am not entirely sure if everything works correctly but the results look promising. I would be very grateful for any confirmation or corrections.