How to address point boundary condition?

1

i want to apply a point boundary condition like the figure show. point load and point support. how to perform it in FEniCS？

image
image1642×921 28.6 KB

2

Hello
Lets say you want to apply a point load equal to 100 N downward in Y direction at coordinate (5,5):

class Load(UserExpression):

    def __init__(self, **kwargs):
        super().__init__(degree=kwargs["degree"])

        self.point = kwargs['point']
        self.value = kwargs['value']

    def eval(self, value, x):

        if near(x[0], self.point[0]) and near(x[1], self.point[1]):
            value[0] = self.value[0]
            value[1] = self.value[1]
        else:  
            value[0] = 0
            value[1] = 0

    def value_shape(self):
        return (2,)

You can call it by:

T = Load(point=(5, 5), value=(0, -100), degree=1)

Regarding applying boundary condition at a specific point:

Lest say you want to fix a point at coordinate (0,0). You can do as following:

def point(x, on_boundary):
    tol = DOLFIN_EPS
    return (abs(x[0] - 0) < tol) and (abs(x[1] -0) < tol)

BC_POINT = DirichletBC(V, Constant((0., 0.)), point,method="pointwise")
1 Like
4

thanks leo, i can figure it out

5

I clumsily made a few changes to extend the ‘Load’ class to a 3d cantilever beam model, and it seemed to work. However, I need to go further and use this functionality to apply the load to multiple points, say, along a boundary, but I’m afraid I don’t know how the class produces the output it does. I’ll explain.

You can see how I modified the above code in the following example:

from __future__ import print_function
from fenics import *
from ufl import nabla_div
import numpy as np

#Updated class definition here
class Load(UserExpression):

    def __init__(self, **kwargs):
        super().__init__(degree=kwargs["degree"])

        self.point = kwargs['point']
        self.value = kwargs['value']

    def eval(self, value, x):

        if near(x[0], self.point[0]) and near(x[1], self.point[1]) and near(x[2], self.point[2]):
            value[0] = self.value[0]
            value[1] = self.value[1]
            value[2] = self.value[2]
        else:  
            value[0] = 0
            value[1] = 0
            value[2] = 0

    def value_shape(self):
        return (3,)

# Scaled variables
l = 2; W = 0.2
mu = 1
rho = 1
delta = W/l
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(l, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Load(point=(2, 0, 0), value=(0, 0, -2), degree=2)
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

When I inspect the variable T, it says:

image
image749×34 1.77 KB

So, three questions:

  1. Where in the class definition or execution of “Load” is the node address accessed in order to create the variable?
  2. I assume that “degree” is the polynomial order of the element. What should inform this choice?
  3. And most importantly, how can I redefine the class to apply the load on multiple nodes on a boundary? In this case, it would be on all nodes along the boundary near x[0]=2 and x[2] = 0.2.

Thanks for your help!

6

Could someone tell me what is the meaning of degree=1 in Load(point=(5, 5), value=(0, -100), degree=1)?
Thank you.