### Problem implementing a point load with Taylor-Hood element

Hello everybody,

I am currently trying to model the problem of an elastic and incompressible plate submitted to a concentrated force on the middle of the top boundary.

To model the concentrated force, I first created and marked a subdomain for the point of application of the force, then I used the point integral dP to implement it in the variational formulation.

Because it is an incompressible material I use the Taylor-Hood element, but it seems to be incompatible with the use of a point integral.

I get indeed the following error:

```
*** Error: Unable to assemble form over vertices.
*** Reason: Expecting test and trial spaces to only have dofs on vertices for point integrals.
*** Where: This error was encountered inside Assembler.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset: 4c59bbdb45b95db2f07f4e3fd8985c098615527f
```

Meaning that I cannot use $P_2\times P_1$`P`_{2}×`P`_{1} element, but $P_1\times P_{...}$`P`_{1}×`P`_{...} . For example, for an element $P_1\times P_1$`P`_{1}×`P`_{1} the compilation runs without problem.

Do you have any solution to make the point integral compatible with Taylor-Hood element?

Else, do you have a better approach to adress this problem?

You will find below a minimal code reproducing the error.

Thank you in advance!

Haythem

```
"""
FEniCS program
The model is used to simulate an elastic material under concentrated load.
"""
from __future__ import division
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
# Mechanical variables
E = 1; nu = 0.3
lambda_ = E*nu/((1+nu)*(1-2*nu))
mu = E/(2*(1+nu))
# Create mesh and define function space
mesh = UnitSquareMesh(50, 50)
# Define function spaces
V = VectorElement('P', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V*Q)
# Define subdomains for boundary conditions
class LeftRight(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and (near(x[0], 0., tol) or near(x[0], 1., tol))
clamped_boundary = LeftRight()
# Define boundary conditions
bc = DirichletBC(W.sub(0), Constant((0, 0)), clamped_boundary)
# Define subdomains for point loads
class PointLoad(SubDomain):
def __init__(self, xcoords):
self.xcoords = xcoords
SubDomain.__init__(self)
def inside(self, x, on_boundary):
tol = 1E-2
return near(x[1], 1., tol) and near(x[0], self.xcoords, tol)
force_markers = VertexFunction("size_t", mesh, 0)
force_markers.set_all(0)
xcoords_loadpoint = 0.5
point_load = PointLoad(xcoords_loadpoint)
point_load.mark(force_markers, 1)
# Define the point integral measure
dP = Measure('dP', domain=mesh, subdomain_data=force_markers)
# Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def epsd(u):
return epsilon(u) - 1./2*nabla_div(u)*Identity(2) # deformations planes
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
d = u.geometric_dimension() # space dimension
pointload_ = Constant((0., -1e-2))
a = (inner(2*mu*epsd(u), epsd(v)) + p*nabla_div(v) + q*nabla_div(u) - 1./lambda_*p*q)*dx
L = dot(pointload_, v)*dP(1)
# Solve variational problem
A = assemble(a)
b = assemble(L)
bc.apply(A, b)
# Compute solution
w = Function(W)
solve(A, w.vector(), b)
# Split the mixed solution using deepcopy
# (needed for further computation on coefficient vector)
(u, p) = w.split(True)
# # Split the mixed solution using a shallow copy
(u, p) = w.split()
# Plot solution
plt.figure()
plot(u, title="velocity")
```