Deformation of a solid for a different force on each node

General
1

Hello everyone, I’m kind of new to fenics.

I´m trying to calculate the deformations of a solid when a force is applied, however I need this force to be different for every node on the boundary to be deformed. Is there a way to define the function f to have a different value for every node?

The forces I must apply look like this.

[[-2.702286116082797e-07, 2.38704010967029e-08],
[-1.5016339088079627e-07, 1.4424024902422665e-07],
[-1.7847667186907969e-06, 1.0988299194180779e-06],
[4.484789160742953e-08, -1.888436087663511e-05],
[-2.604204881798994e-05, -5.9228418740760066e-05],
[1.5268030506832267e-06, -7.803037228975471e-05],
[-2.407123700310234e-06, -7.708842114925285e-05],
[3.034962642204462e-05, -6.165614595515307e-05],
[6.4549847930940455e-06, -1.2834797489998583e-05],
[1.0827009339211713e-06, 5.467043043179808e-06],
[8.337026794717074e-07, -1.8921281666307611e-07],
[1.8334420157521985e-07, 2.990004153563378e-08]]

2

Yes, that is possible.
However, it would be good if you could make a minimal example of your problem.

It is unclear to me how you want to map this 2D array to your mesh. Do you know what coordinate should have what force? Also, are you using dolfinx or legacy dolfin?

3

I have this code that solves my problem for two different domains, I’d like to create one in which I can create n domains, so to speak. Now, I know the coordinates in which the forces must be (they come from another FEM simulation), so I was planning on using the same elements and functions in order for them to fit, so to speak.

Here is the code with two subdomains, using dolfin legacy:

  from dolfin import *
  
  L = 20
  H = 20
  Nx = 250
  Ny = 250
  
  
  mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)
  
  V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
  
  class Sub1(SubDomain):
      def inside(self, x, on_boundary):
          return (between(x[0], (0, 10)) and between(x[1], (0, 1)))
  
  class Sub2(SubDomain):
      def inside(self, x, on_boundary):
          return (between(x[0], (10, 20)) and between(x[1], (0, 1)))
  
  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((1, 1))
  f2 = Constant((0, 0))
  
  dx = Measure("dx")[domains]
  
  
  def eps(v):
      return sym(grad(v))
  
  E = Constant(1e5) #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_streso"
  
  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 left(x, on_boundary):
      return near(x[0], 10)
  
  bc = DirichletBC(V, Constant((0.,0.)), left)
  
  u = Function(V, name="Displacement")
  solve(a == l, u, bc)
  
  plot(1e3*u, mode="displacement")