Applying a force on a face

variational formulation
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

22

Hello, Ihave tried to visualize delta in the code below, could you please try to see it? I am getting an error (output message) in paraview when I open PRO.pvd there and wondering if you will too. Thanks in advance.

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

        def value_shape(self): 
            print("ah")
            return (3, ) 
             

    delta = Delta(eps=1E-4, x0=[-0,-0,10], degree=3)
    PRO = project(delta, V)
    plot(PRO)

    vtkfile_ = File('PRO.pvd')
    vtkfile_ << PRO
23

You still do not supply a full code.
This works for me in paraview 5.7.0:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh, "CG", 1)
class Delta(UserExpression):
   def __init__(self, eps, x0, **kwargs):
      print("init function")
      self.eps = eps
      self.x0 = x0
      UserExpression.__init__(self, **kwargs)
   def eval(self, values, x):

            eps = self.eps
            values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
            values[1] = 0
            #values[2] = -1

   def value_shape(self):
            print("ah")
            return (2, )


delta = Delta(eps=1E-4, x0=[0.2,0.3], degree=3)
PRO = project(delta, V)
plot(PRO)

vtkfile_ = File('PRO.pvd')
vtkfile_ << PRO
24

Thanks a lot, is it possible to convert this example to 3d like this? I am getting this error when I tried I was able to see 3D axis earlier.I don’t know why it did not work now.
NotImplementedError: It is not currently possible to manually set the aspect on 3D axes

from dolfin import *
import numpy as np
mesh = UnitCubeMesh(10,10,10)
V = VectorFunctionSpace(mesh, "CG", 1)
plot(mesh)
class Delta(UserExpression):
   def __init__(self, eps, x0, **kwargs):
      print("init function")
      self.eps = eps
      self.x0 = x0
      UserExpression.__init__(self, **kwargs)
   def eval(self, values, x):

            eps = self.eps
            values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
            values[1] = 0
            values[2] = -1

   def value_shape(self):
            print("ah")
            return (3, )


delta = Delta(eps=1E-4, x0=[0.1,0.5,0.2], degree=3)
PRO = project(delta, V)
plot(PRO)

vtkfile_ = File('PRO.pvd')
vtkfile_ << PRO
25

I will not run your code as you have not formated it properly with ```.
The 3D plotting issue is Fixed in dolfin master, and is because matplotlib changed behavior after the 2019.1.0 release. You can find a working image for plotting with matplotlib on quay.io.
However, in 3D, save it to file and visualize it in Paraview.

26

Ok thank you. I am running the code you provided in jupyter notebook and I can see this plot
image
However, pvd file 1 KB and I am getting error still. What should I do?