Hi. I want to solve the Stokes equations in a unit cube with no slip boundary conditions in all faces of the cube and one point force in the center of the cube (pointing, for example, along the z axis). This is to solve the Green’s function of the Stokes flow. Bellow I’m sharing the code (using legacy dolfin) that I wrote for that purpose but I’m missing the part were I assign a vectorial point force to the center of the cube.

I’ve seen how you can create a PointSource like: delta = PointSource(V, point, 1.0) and then apply it to the rhs like: delta.apply(bb) but that’s not quite what I want to do. I don’t even know what is happening when I apply that point source that is a single scalar. What I want to have is f = (1.0, 0.0, 0.0) but applied to a single point in my domain. A rhs equal to:

f * n_vector * 3D_DiracDelta(x_vector -x0_vector)

where f is a scalar and n_vector is a unit vector representing the orientation of the point force in space.

Here is my code missing the point source term:

```
from dolfin import *
mesh = UnitCubeMesh(10, 10, 10)
# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element()]))
# Boundaries
def all_boundaries(x, on_boundary):
return x[0] > (1.0 - DOLFIN_EPS) or x[0] < (DOLFIN_EPS) or x[1] > (1.0 - DOLFIN_EPS) or x[1] < (DOLFIN_EPS) or x[2] > (1.0 - DOLFIN_EPS) or x[2] < (DOLFIN_EPS)
def top_y(x, on_boundary): return x[1] > (1.0 - DOLFIN_EPS)
noslip = Constant((0.0, 0.0, 0.0))
bc0 = DirichletBC(W.sub(0), noslip, all_boundaries)
bcs = [bc0]
# lid = Constant((1.0, 0.0, 0.0))
# bc1 = DirichletBC(W.sub(0), lid, top_y)
# bcs = [bc0, bc1]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0, 0.0))
a = inner(grad(u), grad(v))*dx + div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx
b = inner(grad(u), grad(v))*dx + p*q*dx
A, bb = assemble_system(a, L, bcs)
P, btmp = assemble_system(b, L, bcs)
solver = KrylovSolver("minres", "amg")
solver.set_operators(A, P)
U = Function(W)
solver.solve(U.vector(), bb)
u, p = U.split()
```

I’m new to Fenics. I started trying to use Fenicsx for this but as a newby everything looked more complicated there and even experienced users were having trouble with “scalar” point sources. Thanks in advance.