Can a subdomain just be a node?

General
1

Hello,

I need to define a subdomain in FeniCs in order for the subdomain to only be one node, as I have a list of forces, each of which I need to apply into one node of a mesh. To simplify, I tried to implement the following code, which should only apply a force on points (0.1,0.) and (0.9,0.9), which by the way the mesh was created, I know are nodes. However, no deformation is occuring. Is there any other way to define this subdomains?

Here is my code:

  mesh = UnitSquareMesh(10,10)
  
  V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
  
  class Sub1(SubDomain):
      def inside(self, x, on_boundary):
          eps = 1e-05
          return (between(x[0], (0.1-eps,0.1+eps)) and between(x[1], (0.1-eps,0.1+eps)) ) 
  
  class Sub2(SubDomain):
      def inside(self, x, on_boundary):
          eps = 1e-05
          return (between(x[0], (0.9-eps,0.9+eps)) and between(x[1], (0.2-eps,0.2+eps)) ) 
  
  sub1 = Sub1()
  sub2 = Sub2()
  
  domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
  domains.set_all(0)
  
  sub1.mark(domains, 1)
  sub2.mark(domains, 2)
  
  f1 = Constant((0, -200000))
  f2 = Constant((0, -10))
  
  dx = Measure("dx")[domains]
  
  
  def eps(v):
      return sym(grad(v))
  
  E = Constant(1) #Young's modulus: Measure of the stiffness of a material and is defined as the ratio of stress to strain in the elastic deformation range.
  nu = Constant(0.3) # measure of the relative contraction or expansion of a material perpendicular to the direction of applied stress
  model = "plane_stres"
  
  mu = E/2/(1+nu)
  lmbda = E*nu/(1+nu)/(1-2*nu)
  if model == "plane_stress":
      lmbda = 2*mu*lmbda/(lmbda+2*mu)
  
  def sigma(v):
      return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
  
  du = TrialFunction(V)
  u_ = TestFunction(V)
  
  a = inner(sigma(du), eps(u_))*dx
  l = inner(f1, u_)*dx(1) + inner(f2, u_)*dx(2)
  
  def up(x, on_boundary):
      return near(x[1], 1)
  
  bc = DirichletBC(V, Constant((0.,0.)), up)
  
  u = Function(V, name="Displacement")
  solve(a == l, u, bc)