Applying a force on a face

variational formulation
1

I am trying to solve a simple 3D elasticity problem, but I can’t apply the force correctly.

I want to fix the displacement in all directions on left_face, and apply a force in -X direction on the right_face.
So I defined the subdomains and want to make the integral of the load f on right_boundary, so I defined ds on that subdomain. But I am not sure if this is correct…

> from __future__ import print_function
> from fenics import *
> from ufl import nabla_div
> import dolfin
> import mshr
> import math
> 
> # creating the mesh
> L = 300
> y = 30
> z = 30
> pt1 = Point(0,0,0)
> pt2 = Point(L,y,z)
> rectangle = mshr.Box(pt1,pt2)
> mesh = mshr.generate_mesh(rectangle,10)
> 
> # defining subdomains
> V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
> 
> class Left_face(SubDomain):
>     def inside(self, x, on_boundary):
>         return on_boundary and near(x[0],0)
>     
> class Right_face(SubDomain):
>     def inside(self, x, on_boundary):
>         return on_boundary and near(x[0],L)
>     
> left_face = Left_face()
> right_face = Right_face()
> 
> sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
> sub_domains.set_all(0)
> 
> left_face.mark(sub_domains, 1)
> right_face.mark(sub_domains, 2)
> 
> ds = ds(subdomain_data = sub_domains, subdomain_id = 2)
> 
> bc = DirichletBC(V, Constant((0.,0.,0.)), bcs = right_face)
> 
> 
> # define structural analysis
> E = Constant(1e5)
> nu = Constant(0.3)
> 
> 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(3) + 2*mu*eps(v)
> 
> f = Constant((-200,0,0)) ##########
> du = TrialFunction(V)
> u_ = TestFunction(V)
> a = inner(sigma(du),eps(u_))*dx
> l = dot(f,u_)*ds(1)  #############
> 
> u = Function(V, name='Displacement')
> solve(a == l, u, bcs=bc)

That error is occurring when I run the code:

> runfile('/home/jiri/Desktop/linux/codigos/untitled0.py', wdir='/home/jiri/Desktop/linux/codigos')
> Traceback (most recent call last):
> 
>   File "<ipython-input-43-121d3ee6abe0>", line 1, in <module>
>     runfile('/home/jiri/Desktop/linux/codigos/untitled0.py', wdir='/home/jiri/Desktop/linux/codigos')
> 
>   File "/usr/lib/python3/dist-packages/spyder/utils/site/sitecustomize.py", line 705, in runfile
>     execfile(filename, namespace)
> 
>   File "/usr/lib/python3/dist-packages/spyder/utils/site/sitecustomize.py", line 102, in execfile
>     exec(compile(f.read(), filename, 'exec'), namespace)
> 
>   File "/home/jiri/Desktop/linux/codigos/untitled0.py", line 48, in <module>
>     bc = DirichletBC(V, Constant((0.,0.,0.)), bcs = right_face)
> 
>   File "/usr/lib/python3/dist-packages/dolfin/fem/dirichletbc.py", line 107, in __init__
>     if isinstance(args[2], cpp.mesh.SubDomain):
> 
> IndexError: tuple index out of range
2

Hi

There’s a mistake in your code:

bc = DirichletBC(V, Constant((0.,0.,0.)), bcs = right_face)

should be

bc = DirichletBC(V, Constant((0.,0.,0.)), sub_domains, 2)
3

Or simply

bc = DirichletBC(V, Constant((0.,0.,0.)), right_face)
4

Hello, I want to have a cantilever beam fixed at one end and apply a point in the vertical direction at the free end in addition to body force , what would be done to do that?

5

You could look at the PointSource class.

6

Thanks, is the BC here applied on a point? I defined subdomains and want to create a cantilever beam
class Left_face(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0],0)

            class Right_face(SubDomain):
                 def inside(self, x, on_boundary):
                        return on_boundary and near(x[0],L)
            left_face = Left_face()
            right_face = Right_face()
            sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
            sub_domains.set_all(0)
             
            left_face.mark(sub_domains, 1)
            right_face.mark(sub_domains, 2)
             
            ds = ds(subdomain_data = sub_domains, subdomain_id = 2)
            bc = DirichletBC(V, Constant((0.,0.,0.)), right_face)

bc is giving an error when I try, itsays that

Error: Unable to create Dirichlet boundary condition.
*** Reason: Expecting a scalar boundary value but given function is vector-valued.
*** Where: This error was encountered inside DirichletBC.cpp.
*** Process: 0

7

As i have told you many times @amy, you Need to supply the a complete code, by that i mean a code of Max 25 lines that reproduces the error, otherwise no one can Give you additional support. The error message indicates that V is a scalar function space, Which you are supplying a 3D vector.

8 
  Ok, Here is the complete code.
    from dolfin import *
    from __future__ import print_function
    from fenics import *
    from ufl import nabla_div
    import dolfin
    import mshr
    import math
    import matplotlib.pyplot as plt
    L = 10.
    W = 1.
    # Mesh for cantilever beam
    mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
    N = 5
    #mesh = UnitSquareMesh(2*N,2*N)
    V = FunctionSpace(mesh,"CG",1)

    # Setting a value using a PointSource, as suggested in the post title:
    p = PointSource(V,Point(1-DOLFIN_EPS,0.5),1.0)
    p.apply(u.vector())

    # Can also set using a BC, which has the advantage of being convenient
    # for applying a constraint during a solve.  The trick is that you need to
    # pass the optional argument "pointwise".
    #bc = DirichletBC(V,Constant(10.0),"near(x[0],0)&&near(x[1],0)&&near(x[2],0)","pointwise")


    # Note: In both cases, the point where the value is set needs to correspond to
    # a DoF of the FunctionSpace.  

    # Check: Function is set to 1 at both points:


    class Left_face(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(x[0],0)
         
    class Right_face(SubDomain):
         def inside(self, x, on_boundary):
                return on_boundary and near(x[0],L)
    left_face = Left_face()
    right_face = Right_face()
    sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    sub_domains.set_all(0)
     
    left_face.mark(sub_domains, 1)
    right_face.mark(sub_domains, 2)
     
    ds = ds(subdomain_data = sub_domains, subdomain_id = 2)
    #bc = DirichletBC(V, Constant((0.,0.,0.)), sub_domains, 2)
    bc = DirichletBC(V, Constant((0.,0.,0.)), right_face)
    #bc = DirichletBC(V,Constant(10.0),"near(x[0],0)&&near(x[1],0)&&near(x[2],0)","pointwise")
    u = Function(V)
    E = Constant(1e5)
    nu = Constant(0.3)
     
    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(3) + 2*mu*eps(v)
    f = Constant((0.0, 0.0, -1e-3))
    #f = Constant((-20,0,0)) ##########
    du = TrialFunction(V)
    u_ = TestFunction(V)
    a = inner(sigma(du),eps(u_))*dx
    l =dot(f, u_)*dx +  dot(u_,p)*ds(1)  #############

    u = Function(V, name='Displacement')
    solve(a == l, u, bc)
9

As @dokken said, when using your FunctionSpace V, the boundary condition has to be scalar, as your error message says:

bc = DirichletBC(V, Constant(0.0), right_face)
1 Like
10

ok thanks Dirichlet means displacement. I didn’t understand how this means a applied force.
Could you please explain it?

from dolfin import *
N = 5
mesh = UnitSquareMesh(2*N,2*N)
V = FunctionSpace(mesh,"CG",1)

# Setting a value using a PointSource, as suggested in the post title:
p = PointSource(V,Point(1-DOLFIN_EPS,0.5),1.0)
p.apply(u.vector())

# Can also set using a BC, which has the advantage of being convenient
# for applying a constraint during a solve.  The trick is that you need to
# pass the optional argument "pointwise".
bc = DirichletBC(V,Constant(1.0),"near(x[0],0)&&near(x[1],0.5)","pointwise")
u = Function(V)
bc.apply(u.vector())

# Note: In both cases, the point where the value is set needs to correspond to
# a DoF of the FunctionSpace.  

# Check: Function is set to 1 at both points:
from matplotlib import pyplot as plt
plot(u)
plt.show()
11

Applied forces at a point has been outlined here: Implement point source in the variational formulation Where you either use a mollified dirac delta or the point source class

12

isn’t it possible using to do it via pointsource then? I did not get an answer for my quaetion in the previous post. Can’t I turn the previous example to 3D case?
Would it be correct if I make it 3D like this
bc = DirichletBC(V,Constant(0.,0.,0.),“near(x[0],1)&&near(x[1],1)&&near(x[2],1)”,“pointwise”)

This gives an error for instance :UFLException: Invalid cell 0.0.

        from dolfin import *
        from dolfin import *
        from __future__ import print_function
        from fenics import *
        from ufl import nabla_div
        import dolfin
        import mshr
        import math
        import matplotlib.pyplot as plt
        N = 5
        L=10
        W=1
        #mesh = UnitSquareMesh(2*N,2*N)
        mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
        V = VectorFunctionSpace(mesh,"CG",1)
        u = Function(V)
        # Setting a value using a PointSource, as suggested in the post title:
        p = PointSource(V,Point(1-DOLFIN_EPS,0.5),1.0)
        p.apply(u.vector())

    # Can also set using a BC, which has the advantage of being convenient
    # for applying a constraint during a solve.  The trick is that you need to
    # pass the optional argument "pointwise".
    bc = DirichletBC(V,Constant(0.,0.,0.),"near(x[0],1)&&near(x[1],1)&&near(x[2],1)","pointwise")

    bc.apply(u.vector())

    # Note: In both cases, the point where the value is set needs to correspond to
    # a DoF of the FunctionSpace.  

    # Check: Function is set to 1 at both points:
    from matplotlib import pyplot as plt
    plot(u)
    plt.show()
14

@amy,
Yet again, you do not supply a working example. additionally, it contains loads of unnecessary code which you have commented out. If you expect people to help you, please make the problems executable and readable for other people.

W.r.t. to your code, you define two functionspaces, U and V. And We cant see the definition of U since you haven’t given a complete code. If U and V are not the same function space, you cannot apply boundary conditions from one space to another.

15

Sorry, there is only one function space, I corrected, this works now but I am not sure if it is correct or not. I wanna have distributed load (body force) along the cantilever and a point force.
from dolfin import *
from ufl import nabla_div

# Load mesh and define function space
L = 10.
W = 1.
# Mesh for cantilever beam
Omega = BoxMesh(Point(0, 0, 0), Point(L, W, W),100, 3, 3)
Omega.init()
V = VectorFunctionSpace(Omega, "CG", 1)
tol = 1.e-8 # tolerance
#Dirichlet on left boundary
def l_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)
class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs) 
    def eval(self, values, x):
        eps = self.eps
        values[0] =-1 
        #eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
        values[1] = 0
        values[2] = 0

    def value_shape(self): return (3, )

delta = Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5)

zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, l_boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, -1e-3))
# Elasticity parameters
E = 1e5
nu = 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Strain
def eps(u):
    return sym(grad(u))
model = "plane_stress"
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(u):
    return lmbda*tr(eps(u))*Identity(3) + 2.0*mu*eps(u)



# Define variational problem

a = inner(sigma(u), eps(v))*dx

l=  dot(f, v)*dx +dot(v,delta)*ds
# Compute solution
u = Function(V)
solve(a == l, u, bc)
print("Maximal deflection:", -u(L,1,1))

I don’t know the functions of values for instance I took that part from here
Implement point source in the variational formulation

16

@dokken Hello,
Can class Delta here work for creating a point force at the tip of the cantilever?

17

@amy you have been supplied with all the necessary information to decide if this works or not. You keep on asking very general questions.
Have you looked at the results when using the Delta function? Does it look reasonable. Save the deformation field to VtK or XDMF, and visualize it in paraview

18

Hello, I have a specific question, in this part of the code I put above, I was thinking x0 is a location for the point force, but when I change the x0, nothing changes in the displacement value. Why would it be?

class Delta(UserExpression):
        def __init__(self, eps, x0, **kwargs):
            self.eps = eps
            self.x0 = x0
            UserExpression.__init__(self, **kwargs) 
        def eval(self, values, x):
            eps = self.eps
            values[0] =-1 
            #eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
            values[1] = 0
            values[2] = 0

        def value_shape(self): return (3, )

    delta = Delta(eps=1E-4, x0=np.array([0.5, 0, 0]), degree=5)
19

To check if the delta function is doing what i should, project it to a CG1 function space, and visualize the result in paraview to see if it changes or not

20

Could you type the code to project it?

21

You Need to learn how to use dolfin and check the source code before you Ask questions.

There are many demos using these normal operations. See for instance this demo for the project function