How to define non-uniform traction boundary conditions at nodes

[French_fries]

Hello everyone，as shown, for the elasticity problem in a 2D area (left side fixed, horizontal traction applied on the right side), if the traction force is constant, the code is as follows:

from fenics import *

mesh = RectangleMesh(Point(0., 0.), Point(5, 5), 10, 10)
V = VectorFunctionSpace(mesh, "P", 1)

def left(x, on_boundary):
    return near(x[0], 0) and on_boundary
boundary_mark = MeshFunction("size_t", mesh, 1)
class Rightboundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and abs(x[0]-5) < tol

right = Rightboundary()
right.mark(boundary_mark, 1)
ds = Measure('ds', domain=mesh, subdomain_data = boundary_mark)

E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

bc = DirichletBC(V, Constant((0, 0)), left)
T = 100
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_)) * dx
l = dot(Constant((T, 0)), u_) * ds(1)

u = Function(V, name="Displacement")
solve(a == l, u, bc)
File("result/U.pvd") << u

If the traction force at each node is different (assuming from bottom to top, they are (T0, 0), (T1, 0), (T2, 0)…), how should it be handled? I hope to receive some help, thank you.

[dokken]

I assume that the traction should be continuous between the nodes, maybe linear? If so, create a first order Lagrange space, and assign the traction to the relevant node, using V.tabulate_dof_coordinates() and a distance computation from your input points as shown in

How to represent non-uniform first and second type boundary conditions?

First of all, your coordinates make no sense, as will never match with Please pay special attention to details before posting code. Your code is also way to complicated, introducing boundary meshes, sub meshes of boundary meshes etc. Here is a minimal working example of what you want to achieve: from dolfin import * import numpy as np mesh = UnitSquareMesh(4, 4) VT = FunctionSpace(mesh, "CG", 1) T_data = np.array([0, 50, 100, 150, 0]) T_coordinates = np.array([[0, 0], [0., 0.25], [0., 0.…

[French_fries]

In fact, I’ve already made a similar attempt before:

from dolfin import *
import numpy as np

# Create mesh and function space
mesh = RectangleMesh(Point(0., 0.), Point(5, 5), 5, 5)
V = VectorFunctionSpace(mesh, "CG", 1)

T_data = [0, 50, 100, 120, 150, 0]
T_coordinates = [[5, 0], [5, 1], [5, 2], [5, 3], [5, 4], [5, 5]]

class Rightboundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 5) and on_boundary

right = Rightboundary()
vertex_mark = MeshFunction("size_t", mesh, 0, 0)
right.mark(vertex_mark, 1)
vertices = vertex_mark.where_equal(1)
vertex_to_dof_map = vertex_to_dof_map(V)

T_load = Function(V)

# Get all boundary coordinates
dof_coords = V.tabulate_dof_coordinates()
boundary_dofs = [vertex_to_dof_map[vertex] for vertex in vertices]
boundary_coords = np.array([dof_coords[boundary_dof] for boundary_dof in boundary_dofs])

points_A = np.expand_dims(boundary_coords, 1)
points_B = np.expand_dims(T_coordinates, 0)
distances = np.sum(np.square(points_A - points_B), axis=2)
is_close = distances < 1e2 * DOLFIN_EPS
positions = np.nonzero(is_close)
for row, col in zip(*positions):
    T_load.vector()[boundary_dofs[row]] = T_data[col]
    print(row, col, boundary_coords[row], T_coordinates[col], T_data[col])

E = Constant(50e3)
nu = Constant(0.2)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

boundary_mark = MeshFunction("size_t", mesh,1)
right.mark(boundary_mark, 1)
ds = Measure('ds', domain=mesh, subdomain_data = boundary_mark)

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_)) * dx
l = dot(T_load, u_) * ds(1)

u = Function(V, name="Displacement")
solve(a == l, u)

However, the output results were missing some parts：

1 0 [5. 0.] [5, 0] 0
3 1 [5. 1.] [5, 1] 50
5 2 [5. 2.] [5, 2] 100
Solving linear variational problem.

On the other hand, as far as I know, traction should be defined in the form of a 2D vector because uniform constant traction is defined as T =Constant((0，0)). After I modified the form of the traction:
T_data = [[0,0], [50,0], [100,0], [120,0], [150,0], [0,0]]
the following error occurred:

Traceback (most recent call last):
  File "/home/dyfluid/Desktop/contact2D/node.py", line 34, in <module>
    T_load.vector()[boundary_dofs[row]] = T_data[col]
    ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
IndexError: Indices to set must be a 1D array

I hope you can provide some help, thank you.

[dokken]

@French_fries

T_load = Function(V)

Make T_load in a scalar space, i.e.
T_space = FunctionSpace(mesh, "Lagrange", 1),
T_load = Function(T_space), then in turn assign the T_data.

Then, instead of T_load in the variational form, use:

@French_fries

l = dot(T_load, u_) * ds(1)

Traction = ufl.as_vector((T_load,0))
l = dot(Traction, u_)*ds(1)

[French_fries]

Thank you for your guidance, my issue has been resolved. I have a question regarding the definition of traction: in the FEniCS tutorial, T is defined as the force along the boundary normal. So, for a vertical boundary (x=0), does the traction condition have to be set horizontally (T = Constant(a, 0))? If it is set as T = Constant(a, b), what would be the result? I hope you can answer my question, thank you.

[dokken]

\sigma \cdot n should be zero in the tangential direction if you are only applying normal stresses.
See for instance page 3 of:

elasticity.pdf

In particular:

Stress at a point x can be mathematically defined as limA→{x} tA/A, where tA repre-
sents the internal force applied across the cutting surface A, which encloses x and shrinks
to the point x. The force tA is a vector and does not necessarily have to be orthogonal to
the surface. The component of this force that is normal to the surface is known as the nor-
mal stress (σn), and the component tangential to the surface is termed the shearing stress
(τn), with the subscript n indicating a normal vector to the cutting surface.

and Page 6:

2. Neumann boundary conditions (also known as traction boundary conditions).
   These conditions specify the traction vector t (force per unit area) on the Neumann bound-
   ary ΓN . Mathematically: σ · n = t on ΓN , where σ is the stress tensor, n is the unit
   outward normal vector to the boundary, and t is a given function. Physically, this means
   that the boundaryΓN is subjected to a specified force or traction

[French_fries]

Sure, I understand now. Thank you for your patient explanation!