UPDATED: Constructing nodal force vector involving shape function derivatives

2

[POST UPDATED FOR CLARITY, SHORTER MWE]

I construct a vector of nodal forces, which are functions of the coordinates of some points, given as a list. If a cell contains one of these points, it is to be assigned nodal forces. To simplify this post, without a loss of generality, these nodal dof forces are the element shape function derivatives (“strains”), evaluated at the point’s coordinates. Example code follows. I am looking for a simpler way to do this, preferably at a higher level using ufl directly (the ideal solution would play nice with dolfin-adjoint “out of the box”).

(For more context: Low level basis function operations in UFL)

from dolfin import *
import numpy as np
from ufl import JacobianInverse

mesh = UnitSquareMesh(1,1)
tree = mesh.bounding_box_tree()
V = FunctionSpace(mesh, 'CG', 1)

# Evaluate the i-th basis function n-th derivative in direction j at point x:
def basis_func_derivs(i, j, x, n=1):
    cell_index = tree.compute_first_entity_collision(Point(*x))
    cell_global_dofs = V.dofmap().cell_dofs(cell_index)
    for local_dof in range(len(cell_global_dofs)):
        if(i==cell_global_dofs[local_dof]):
            cell = Cell(mesh, cell_index)
            values = np.array([0,])

            return V.element().evaluate_basis_derivatives(local_dof, n, x,
                                              cell.get_vertex_coordinates(),
                                              cell.orientation())[j]

    # If none of this cell's shape functions map to the i-th basis function,
    # then the i-th basis function is zero at x.
    return 0.0

# Build a Python list of basis function derivatives dN/dxi:
list_of_basis_funcs_derivs = []
for i in range(0,V.dim()):
    list_of_basis_funcs_derivs += [(lambda i: lambda x : basis_func_derivs(i,0,x,n=1))(i),]

J_inv = project(JacobianInverse(mesh))

def f_i(i,x):      # the nodal force term
    return list_of_basis_funcs_derivs[i](x)*J_inv(x).reshape(2,2)[0,0]

# Some points to evaluate the source term:
pts = [[1,0],[0,1]]
f_vec = np.zeros(V.dim())

for pt in pts:
    f_vec += np.array([f_i(i,pt) for i in range(V.dim())])

print(f_vec)
f = Function(V)
f.vector()[:] = f_vec