Point load and body force

variational formulation
1 
from dolfin import *
import mshr
import numpy as np   
from __future__ import print_function
from fenics import *
from __future__ import print_function
from dolfin import *
from __future__ import division
import matplotlib.pyplot as plt
from ufl import nabla_div
from mshr import *
from  dolfin import mesh as module
# Load mesh and define function space
L = 8.
W = 1.
n0 = int(L*6)
n1 = int(W*6)
n2 = int(W*6)
# Mesh for cantilever beam
Omega = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
# = Box(0., 0., 0., d0, d1, d2, n0, n1, n2)
Omega.init()
V = VectorFunctionSpace(Omega, "CG", 1)
tol = 1.e-8 # tolerance 
# Left boundary 
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
# Right boundary
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 8.0)      
# Point (8,0.5,0.5) at which load is added
class point(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],8.) and near(x[1],0.5,1e-2) and near(x[2],0.5,1e-2)
 #Boundary segments
left = Left()
pt = point()
right = Right()
boundaries = MeshFunction("size_t", Omega,0)
boundaries.set_all(0)
left.mark(boundaries, 1)
pt.mark(boundaries, 2) 
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, 00.0))
# Elasticity parameters
E = 210000.0
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(v))
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(2) + 2.0*mu*eps(u)
ds = ds(subdomain_data=boundaries)
# Define variational problem
ds= Measure("ds")(subdomain_data=boundaries)
g_T= Constant((0.0, 0.0, -0.3))
a = inner(sigma(u), eps(v))*dx
left_end_displacement = Constant((0.0, 0.0, 0.0))
#bc = [DirichletBC(V, left_end_displacement, left)] 
bc = [DirichletBC(V, left_end_displacement, left),DirichletBC(V, g_T, pt, method="pointwise")] 
L=  dot(f, v)*dx +dot(v,g_T)*ds(2)
# Compute solution
u = Function(V)
solve(a == L, u, bc);`
2

hello can anyone help me to fix
this code, there is an error after running a,variational form is wrong I guess

3

As far as i can see, you are doing a 3D problem, and then you have to use Identity(3) in the definition of sigma.

4

I am new on FEM and fenics and modified an example I found on the internet. I am still getting an error, I guess there is another mistake in the definition of l

5

It would help if you post the actual error message

6

RuntimeError Traceback (most recent call last)
in ()
1 u = Function(V)
----> 2 solve(a ==l, u, bc);

/usr/lib/python3/dist-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
218 # tolerance)
219 elif isinstance(args[0], ufl.classes.Equation):
–> 220 _solve_varproblem(*args, **kwargs)
221
222 # Default case, just call the wrapped C++ solve function

/usr/lib/python3/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
245 solver = LinearVariationalSolver(problem)
246 solver.parameters.update(solver_parameters)
–> 247 solver.solve()
248
249 # Solve nonlinear variational problem

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** fenics-support@googlegroups.com

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to extract mesh from form.
*** Reason: Non-matching meshes for function spaces and/or measures.
*** Where: This error was encountered inside Form.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

7

I am running this code now, and getting the error above

from dolfin import *
import mshr
import numpy as np
from future import print_function
from fenics import *
from future import print_function
from dolfin import *
from future import division
import matplotlib.pyplot as plt
from ufl import nabla_div

from mshr import *
from dolfin import mesh as module

Load mesh and define function space

L = 8.
W = 1.

n0 = int(L6)
n1 = int(W6)
n2 = int(W*6)

Mesh for cantilever beam

Omega = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)

= Box(0., 0., 0., d0, d1, d2, n0, n1, n2)

Omega.init()

V = VectorFunctionSpace(Omega, “CG”, 1)

tol = 1.e-8 # tolerance

Left boundary

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

Right boundary

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 8.0)

Point (8,0.5,0.5) at which load is added

class point(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],8.) and near(x[1],0.5,1e-2) and near(x[2],0.5,1e-2)

#Boundary segments
left = Left()
pt = point()
right = Right()

Initialize mesh function for boundary domains

#boundaries = MeshFunction(“size_t”, Omega, file)#“size_t”, mesh, 2, mesh.domains()
#boundaries.set_all(0)
#left.mark(boundaries, 1)
#pt.mark(boundaries, 2)
boundaries = MeshFunction(“size_t”, Omega,0) #mesh omega yerinde kalmisti
#boundaries = CellFunction(‘size_t’, mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
pt.mark(boundaries, 2)
#all_boundary.mark(boundaries, 1)

Define variational problem

f = Constant((0.0, 0.0, 00.0))
E = 210000.0
nu = 0.3

mu, lmbda = Constant(E/(2*(1 + nu))), Constant(Enu/((1 + nu)(1 - 2*nu)))

Strain

def eps(v):
return sym(grad(v))
model = “plane_stress”

mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)
if model == “plane_stress”:
lmbda = 2mulmbda/(lmbda+2*mu)

def sigma(v):
return lmbdatr(eps(v))Identity(3) + 2.0mueps(v)

ds = ds(subdomain_data=boundaries)

Define variational problem

ds= Measure(“ds”)(subdomain_data=boundaries)
g_T= Constant((0.0, 0.0, -0.3))
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
left_end_displacement = Constant((0.0, 0.0, 0.0))

#bc = [DirichletBC(V, left_end_displacement, left)]
bc = [DirichletBC(V, left_end_displacement, left),DirichletBC(V, g_T, pt, method=“pointwise”)]
l= dot(f, v)*dx +dot(v,g_T)*ds(2)
u = Function(V)
solve(a ==l, u, bc);

8

Your code was riddled with typos. Please be more careful when writing your code. The first error is because you did not supply the domain when defining the measure. The next error was due to the fact that you use v inside def eps, which should be a function of u.
The following cleaned up code runs:

from dolfin import *
from ufl import nabla_div

# Load mesh and define function space
L = 8.
W = 1.
# Mesh for cantilever beam
Omega = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
Omega.init()
V = VectorFunctionSpace(Omega, "CG", 1)
tol = 1.e-8 # tolerance
# Left boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
# Right boundary
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 8.0)
# Point (8,0.5,0.5) at which load is added
class point(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],8.) and near(x[1],0.5,1e-2) and near(x[2],0.5,1e-2)
 #Boundary segments
left = Left()
pt = point()
right = Right()
boundaries = MeshFunction("size_t", Omega,0)
boundaries.set_all(0)
left.mark(boundaries, 1)
pt.mark(boundaries, 2)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, 00.0))
# Elasticity parameters
E = 210000.0
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

ds= Measure("ds", domain=Omega, subdomain_data=boundaries)
g_T= Constant((0.0, 0.0, -0.3))
a = inner(sigma(u), grad(v))*dx

left_end_displacement = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, left),DirichletBC(V, g_T, pt, method="pointwise")]
L=  dot(f, v)*dx +dot(v,g_T)*ds(2)
# Compute solution
u = Function(V)
solve(a == L, u, bc)
9

Thank you so much, and how should I learn the displacement at the free end also plot this deformation

10

There are many demos at bitbucket. For instance demo for elasticity.
Instead of using matplotlib you can also save the function and visualize it in Paraview

11

I will look at it thank you

12

Would you tell me ,
-u(10,1,1)[2] is this this command correct? What does [2] mean here ?

13

It’s the z component of the deformation

14

which is the displacement then, right

15

My other question is: In this code if I put an end load, I am not getting logical results. My deflection formula is the superposition of point load and distributed load. I am expecting better prediction of displacement when I increase the mesh.
δ=(qL^4)/8EI+(FL^3)/3EIδ=(qL^4)/8EI+(FL^3)/3EI
Where am I doing a mistake?

16

You need to make a minimal code that illustrates how the behavior is different from what you expect.

17

actually, I am questioning this superposition formula, it is true for Euler Beam right

18

Here for instance I am not considering the body force, and FL^3/3EI gives 0,032.
This code gives 0.79 with 100 elements and with 10 elements it is 0.17

from dolfin import *
from ufl import nabla_div

Load mesh and define function space

L = 20.
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

Left boundary

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

Right boundary

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], L)

Point (8,0.5,0.5) at which load is added

class point(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],L) and near(x[1],0.5,1e-2) and near(x[2],0.5,1e-2)
#Boundary segments
left = Left()
pt = point()
right = Right()
boundaries = MeshFunction(“size_t”, Omega,0)
boundaries.set_all(0)
left.mark(boundaries, 1)
pt.mark(boundaries, 2)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0.0, -0))

Elasticity parameters

E = 1e5
nu = 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(Enu/((1 + nu)(1 - 2*nu)))

Strain

def eps(u):
return sym(grad(u))
model = “plane_stress”
mu = E/2/(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)
if model == “plane_stress”:
lmbda = 2mulmbda/(lmbda+2mu)
def sigma(u):
return lmbdatr(eps(u))Identity(3) + 2.0mu*eps(u)

Define variational problem

ds= Measure(“ds”, domain=Omega, subdomain_data=boundaries)
g_T= Constant((0.0, 0.0, -1e-1))
a = inner(sigma(u), grad(v))*dx

left_end_displacement = Constant((0.0, 0.0, 0.0))
bc = [DirichletBC(V, left_end_displacement, left),DirichletBC(V, g_T, pt, method=“pointwise”)]
l= dot(f, v)*dx +dot(v,g_T)*ds

Compute solution

u = Function(V)
solve(a == l, u, bc)
print(“Maximal deflection:”, -u(20,1,1)[2])

19

I am tring to use quadratic elements instead of linear elements. if I change 1 to 2 in this code, would it be correct
V = VectorFunctionSpace(Omega, “CG”, 1)
V = VectorFunctionSpace(Omega, “CG”, 2)

20

Yes. This is a trivial question. Please read through the FEniCS tutorial which will guide you