I am getting started to FEniCS, and I have a few time to implement my method.
I want to do a structural analysis (fem) of a mesh that I am generating with mshr and dolfin.
It’s a 3D lattice structure (photo attached), the boundary condition is [displacement in y==0] = 0 and it has a punctual force in -y applied in it’s top (x==0,y==Y,z==0), where Y is the size of the structure.
I don’t know how can I deal with that problem. I tried to look some examples, but it didn’t work.
Please, help me!
So i guess you would like to solve a linear elasticity equation with a pointwise source function. You can create thatas described in Pointsource on specific points python.
On the no displacement boundary, you Need to apply a traditional Dirichlet condition on the sub space as described here Rigid Plate Boundary Condition
I tried the code below, but it didn’t work.
I think the problem is in the definition of the boundary conditions.
mu = 1 beta = 1.25 lambda_ = beta force = -2 # Newtons # define space function V = VectorFunctionSpace(mesh,'P', 1) # boundary conditions tol = 1E-14 def clamped_boundary(x, on_boundary): return on_boundary and x[1] < tol bc = DirichletBC(V, Constant(0.), clamped_boundary) # deine 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(3) + 2*mu*epsilon(u) # define more variables u = TrialFunction(V) d = u.geometric_dimension() v = TestFunction(V) f = Constant((0,force,0)) T = Constant((0,0,0)) a = inner(sigma(u),epsilon(v))*dx L = dot(f,v)*dx + dot(T,v)*ds # solve the problem u = Function(V) solve(a == L, u, bc) # print the solution plot(u, title='Displacement', mode='displacement') # print stress s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) von_Mises = sqrt(3./2*inner(s,s)) V = FunctionSpace(malha,'P',1) von_Mises = project(von_Mises, V) plot(von_Mises, title='stress intensity') # strain computing u_magnitude = sqrt(dot(u,u)) u_magnitude = project(u_magnitude, V) plot(u_magnitude, 'Displacement magnitude') print('min/max u:', u_magnitude.vector().min(), u_magnitude.vector().max()) # save the files File('elasticity/displacement.pvd') << u File('elasticity/von_mises.pvd') << von_Mises File('elasticity/magnitude.pvd') << u_magnitude
The problem is that your source f is not a point source, but a constant everywhere.
If you follow Pointsource on specific points python, you should do the following (I’ve set Y=1):
bc_point = DirichletBC(V,Constant(1.0),"near(x[0],0)&&near(x[1],1)&&near(x[2],0)","pointwise")
f = Function(V)
bc_point.apply(f.vector())
You can plot f to now verify that it is indeed a point source.
I think the force is ok now. But this error is still happening:
File "<ipython-input-7-2ac9940c1238>", line 18, in <module> bc = DirichletBC(V, Constant(0.), clamped_boundary) File "/usr/lib/python3/dist-packages/dolfin/fem/dirichletbc.py", line 131, in __init__ super().__init__(*args) RuntimeError:
I don’t know why.
Have you tried
bc = DirichletBC(V, Constant((0.0,0.0,0.0)), clamped_boundary)
I tried it now, but a different error happened…
> File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 220, in solve
> _solve_varproblem(*args, **kwargs)
>
> File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 247, in _solve_varproblem
> solver.solve()
>
> RuntimeError:Post a MWE and the complete error message.
As @volkerk says, please post a Minimal working example and the full error message.
However, I think your boundary condition should be:
bc = DirichletBC(V.sub(1), Constant(0), clamped_boundary)
Is this enforces a zero deformation condition only in the y-direction (sub-space 1).
I tried to remove the force and apply 2 boundary conditions, which the second one (top) is to simulate the effect of a load.
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,0.0))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du),eps(u_))*dx
l = inner(f, u_)*dx
def basebc(x, on_boundary):
return near(x[1],0.) and on_boundary
desloc_top = -0.1
def topbc(x,on_boundary):
return near(x[1],Y) and on_boundary
bc_base = DirichletBC(V.sub(1), Constant(0.), basebc)
bc_top = DirichletBC(V.sub(1), Constant(Y+desloc_top), topbc)
u = Function(V, name='Displacement')
solve(a == l, u, bcs=[bc_base, bc_top])
The idea is that in y == Y the displacement is -0.1.
That’s the error message:
Traceback (most recent call last):
File "<ipython-input-5-a541056cbc14>", line 31, in <module>
solve(a == l, u, bcs=[bc_base, bc_top])
File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 220, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/python3/dist-packages/dolfin/fem/solving.py", line 247, in _solve_varproblem
solver.solve()
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 successfully call PETSc function 'MatXIJSetPreallocation'.
*** Reason: PETSc error code is: 1 ((null)).
*** Where: This error was encountered inside /build/dolfin-LDAJTl/dolfin-2019.1.0/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------As you have not supplied to mesh, I’ve changed the problem to a 2D problem.
If you would like to to the point force (or any pointwise boundary condition), you need to remove the on_boundary argument from your topbc function, and add pointwise to the topbc. I’ve illustrated how to use a pointforce in the code below:
from dolfin import *
mesh = UnitSquareMesh(10,10)
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(2) + 2*mu*eps(v)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
f = Function(V)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du),eps(u_))*dx
l = inner(f, u_)*dx
def basebc(x, on_boundary):
return near(x[1],0.) and on_boundary
X = 0.5
desloc_top = -1e6
def topbc(x,on_boundary):
return near(x[1],1) and near(x[0], X)
bc_base = DirichletBC(V.sub(1), Constant(0.), basebc)
bc_top = DirichletBC(V.sub(1), Constant(desloc_top), topbc,"pointwise")
bc_top.apply(f.vector())
File("f.pvd") << f
u = Function(V, name='Displacement')
solve(a == l, u, bcs=[bc_base])
File("uh.pvd") << u
To use a point deformation, change the size of desloc_top to something reasonable (like 0.05), and add the bc back into the solve command:
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
f = Constant((0,0))
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du),eps(u_))*dx
l = inner(f, u_)*dx
def basebc(x, on_boundary):
return near(x[1],0.) and on_boundary
X = 0.5
desloc_top = -0.05
def topbc(x,on_boundary):
return near(x[1],1) and near(x[0], X)
bc_base = DirichletBC(V.sub(1), Constant(0.), basebc)
bc_top = DirichletBC(V.sub(1), Constant(desloc_top), topbc,"pointwise")
u = Function(V, name='Displacement')
solve(a == l, u, bcs=[bc_base, bc_top])
File("uh.pvd") << u
Sorry, I forgot to show my mesh.
It has around 772533 cells, so think that’s why the solver not working.
It’s the same mesh I attached a photo before
> # size of the mesh
> X = 2
> # Y and Z are calculated from X and nx,ny,nz
>
> # number of cells in each direction
> nx = 2
> ny = 2
> nz = 2
>
> # angle between the cylinders (degrees)
> theta_y = 45
> theta_z = 45
> theta_y = math.radians(theta_y)
> theta_z = math.radians(theta_z)
>
> # distance between two cells
> dxis = X/(nx - 1)
> dy = dxis*math.tan(theta_y)/2
> dz = dxis*math.tan(theta_z)/2
>
> Y = 2*ny*dy
> Z = 2*nz*dz
> # radius of the cylinders
> t = 0.25
>
> # quality of the mesh
> n_malhas = 40
>
> # iniciate usefull variables
> pi = math.pi
>
> # create the standard structure
> base = Point(0,0,0)
>
> pt1_y = Point(dxis/2,dy,0)
> pt2_y = Point(-dxis/2,dy,0)
> pt1_z = Point(0,dy,dz)
> pt2_z = Point(0,dy,-dz)
> pt2_x = Point(0,2*dy,0)
>
> cili1_xy = mshr.Cylinder(base,pt1_y,t,t)
> cili2_xy = mshr.Cylinder(base,pt2_y,t,t)
> cili3_xy = mshr.Cylinder(pt1_y,pt2_x,t,t)
> cili4_xy = mshr.Cylinder(pt2_y,pt2_x,t,t)
> cili1_xz = mshr.Cylinder(base,pt1_z,t,t)
> cili2_xz = mshr.Cylinder(base,pt2_z,t,t)
> cili3_xz = mshr.Cylinder(pt1_z,pt2_x,t,t)
> cili4_xz = mshr.Cylinder(pt2_z,pt2_x,t,t)
>
> cili1_yz = mshr.Cylinder(pt1_y,pt1_z,t,t)
> cili2_yz = mshr.Cylinder(pt1_y,pt2_z,t,t)
> cili3_yz = mshr.Cylinder(pt2_y,pt1_z,t,t)
> cili4_yz = mshr.Cylinder(pt2_y,pt2_z,t,t)
> esfera1 = mshr.Sphere(base,t)
> esfera2 = mshr.Sphere(pt1_y,t)
> esfera3 = mshr.Sphere(pt2_y,t)
> esfera4 = mshr.Sphere(pt1_z,t)
> esfera5 = mshr.Sphere(pt2_z,t)
> esfera6 = mshr.Sphere(pt2_x,t)
> est_padrao = cili1_xy + cili2_xy + cili3_xy + cili4_xy
> est_padrao += cili1_xz + cili2_xz + cili3_xz + cili4_xz
> est_padrao += cili1_yz + cili2_yz + cili3_yz + cili4_yz
> est_padrao += esfera1 + esfera2 + esfera3 + esfera4 + esfera5 + esfera6
> objeto = est_padrao
>
> # loop in X direction
> for i in range(0,nx):
> x = i*dxis
>
> pt = Point(x,0,0)
>
> est = mshr.CSGTranslation(est_padrao,pt)
>
> if i == 0:
> pass
> else:
> objeto += est
>
> est_padrao = objeto
>
> # loop in Y direction
> for j in range(0,ny):
> y = j*2*dy
>
> pt = Point(0,y,0)
>
> est = mshr.CSGTranslation(est_padrao,pt)
>
> if j == 0:
> pass
> else:
> objeto += est
>
> est_padrao = objeto
>
> # loop in Z direction
> for k in range(0,nz):
> z = k*2*dz
>
> pt = Point(0,0,z)
>
> est = mshr.CSGTranslation(est_padrao,pt)
>
> if k == 0:
> pass
> else:
> objeto += est
>
> # generate the mesh
> mesh = mshr.generate_mesh(objeto,n_malhas,'tetgen')
And another question: is there any argument instead of ‘pointwise’ that apply the boundary condition to a plan(XZ in Y = 1, for example)?
Hi, try to make a 2D version of your problem to ensure you are solving the problem in the correct fashion before doing a huge computation. To apply bcs on a full plane, you do not Need an additional argument, as you are applying the condition to more than one point.
Hello I have implemented your suggestions in 3D cantilever to apply pointforce at the free end in adition to body force but the displacement value I get is not like I expected. There must be a mistake how can I fix this?
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
# Left boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary
# Right boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 10.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],10.) and near(x[1],0.5,1e-2) and near(x[2],1-tol,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, -1.0))
# 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
ds= Measure("ds", domain=Omega, subdomain_data=boundaries)
g_T= Constant((0.0, 0.0, -1))
a = inner(sigma(u), eps(v))*dx
left_end_displacement = Constant((0.0, 0.0, 0.0))
bc_base = DirichletBC(V, left_end_displacement, left)
bc_p = DirichletBC(V, g_T, pt, method="pointwise")
L= dot(f, v)*dx +dot(v,g_T)*ds
# Compute solution
u = Function(V)
bc_p.apply(u.vector())
solve(a == L, u, bc_base)
print("Maximal deflection:", -u(10,0.5,0.5)[2])You should apply bc_p to the variational problem, as you do with bc_base.
Thank you, I have tried this commnet below it but there is no change in displacement.
solve(a == L, u, bcs=[bc_base,bc_p])
I had implemented your example above in January 31.For the one with pointforce, you did not provide bc_top = DirichletBC(V.sub(1), Constant(desloc_top), topbc,“pointwise”) to the solve comment as well though.
This is because you were mixing in tolerance inside your subdomain definition, using near and 1-tol. The following code works.
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, 4, 4)
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) and on_boundary
# Right boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 10.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],10.) and near(x[1],0.5,1e-2) and near(x[2],1,1e-2)
#Boundary segments
left = Left()
pt = point()
right = Right()
boundaries = MeshFunction("size_t", Omega, Omega.topology().dim()-1,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, -1.0))
# 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
ds= Measure("ds", domain=Omega, subdomain_data=boundaries)
g_T= Constant((0, 0, -1))
a = inner(sigma(u), eps(v))*dx
left_end_displacement = Constant((0.0, 0.0, 0.0))
bc_base = DirichletBC(V, left_end_displacement, left)
bc_p = 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_p,bc_base])
print("Maximal deflection:", -u(10,0.5,0.5)[2])
If I want to have a point deformation rather than a point force in the end of the cantilever, what should I change?
Right now you have a point deformation. Point force is Added with dirac delta functions.
In the post in January 31, there are two example you had shown,one is for point force one is for point deformations as far as I understand. Then that example is not correct. Is dirac delta the only way to do it?