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?
For instance Here I have delta for point force, based on my understanding value is the force x0 is the position of force applied. I am expecting a change in the displacement when I change x0 but there is no change. Is my understanding is correct? if not how can I fix this code?
from dolfin import *
from ufl import nabla_div
import numpy as np
from dolfin import *
# 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):
print("init function")
self.eps = eps
self.x0 = x0
UserExpression.__init__(self, **kwargs)
def eval(self, values, x):
eps = self.eps
values[0] = 0
values[1] = 0
values[2] = -1
def value_shape(self):
return (3, )
delta = Delta(eps=1E-4, x0=[-0,-0.,-10], degree=3)
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., -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
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,0.5,0.5)[2])A Dirac delta function is a smoothed point load. If you do as I said on January 31st you will get a point load as well.
I the example below, you are not using a delta function
You are setting the function to constant -1 everywhere. Please read through your code carefully.
thank you, it would be great if you can show how to do, I could not find how to fix it
I am gonna stop replying to you, as you are clearly not putting in any effort to solve this yourself. I have already shown you how to do this in Applying a force on a face
Yes, there is this command there for value, for 2d case :
eps/pi/(np.linalg.norm(x-self.x0)2 + eps2)
from dolfin import *
from ufl import nabla_div
import numpy as np
from dolfin import *
# 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):
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] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
def value_shape(self):
print("ah")
return (3, )
delta = Delta(eps=1E-4, x0=[-10,-0.,-0], degree=3)
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., -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
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,0.5,0.5)[2])
but still displacement not affected much by the point (x0 )
Increase the magnitude of the dirac delta function then. A point source in 3D is less effective than 2D
if I want to apply a specific value for the force is there a way to do it? what is the the magnitude of force here