# Help a starter!

**URL:** <https://fenicsproject.discourse.group/t/help-a-starter/2324>\
**Category:** variational formulation\
**Created:** [January 30, 2020, 10:02am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324 "2020-01-30T10:02:52Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Thiago](https://avatars.discourse-cdn.com/v4/letter/t/ec9cab/32.png) [@Thiago](https://fenicsproject.discourse.group/u/Thiago)\
**Post date:** [January 30, 2020, 10:02am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/1 "2020-01-30T10:02:52Z")

</div>

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.

 ![2020-01-29-140058_1920x1080_scrot](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/dbf38d902532a8d894793972bd4e43bea5624782.png)

Please, help me!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 30, 2020, 10:14am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/2 "2020-01-30T10:14:03Z")

</div>

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](https://fenicsproject.discourse.group/t/pointsource-on-specific-points-python/1967).  
On the no displacement boundary, you Need to apply a traditional Dirichlet condition on the sub space as described here [Rigid Plate Boundary Condition](https://fenicsproject.discourse.group/t/rigid-plate-boundary-condition/2056/2)

---

<div class="post-metadata">

**Author:** ![Thiago](https://avatars.discourse-cdn.com/v4/letter/t/ec9cab/32.png) [@Thiago](https://fenicsproject.discourse.group/u/Thiago)\
**Post date:** [January 30, 2020, 10:31am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/3 "2020-01-30T10:31:29Z")

</div>

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
> 
> ```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 30, 2020, 10:40am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/4 "2020-01-30T10:40:41Z")

</div>

The problem is that your source `f` is not a point source, but a constant everywhere.  
If you follow [Pointsource on specific points python](https://fenicsproject.discourse.group/t/pointsource-on-specific-points-python/1967), you should do the following (I’ve set `Y=1`):

```auto
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.

---

<div class="post-metadata">

**Author:** ![Thiago](https://avatars.discourse-cdn.com/v4/letter/t/ec9cab/32.png) [@Thiago](https://fenicsproject.discourse.group/u/Thiago)\
**Post date:** [January 30, 2020, 12:06pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/5 "2020-01-30T12:06:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![volkerk](https://avatars.discourse-cdn.com/v4/letter/v/9d8465/32.png) [@volkerk](https://fenicsproject.discourse.group/u/volkerk)\
**Post date:** [January 30, 2020, 12:34pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/6 "2020-01-30T12:34:17Z")

</div>

Have you tried

```python
bc = DirichletBC(V, Constant((0.0,0.0,0.0)), clamped_boundary)

```

---

<div class="post-metadata">

**Author:** ![Thiago](https://avatars.discourse-cdn.com/v4/letter/t/ec9cab/32.png) [@Thiago](https://fenicsproject.discourse.group/u/Thiago)\
**Post date:** [January 30, 2020, 1:32pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/7 "2020-01-30T13:32:56Z")

</div>

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:
```

---

<div class="post-metadata">

**Author:** ![volkerk](https://avatars.discourse-cdn.com/v4/letter/v/9d8465/32.png) [@volkerk](https://fenicsproject.discourse.group/u/volkerk)\
**Post date:** [January 30, 2020, 1:46pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/8 "2020-01-30T13:46:38Z")

</div>

Post a MWE and the complete error message.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 30, 2020, 1:54pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/9 "2020-01-30T13:54:45Z")

</div>

> [@Thiago](#):
>
> ```auto
> 
> ```

As @volkerk says, please post a Minimal working example and the full error message.  
However, I think your boundary condition should be:

```auto
bc = DirichletBC(V.sub(1), Constant(0), clamped_boundary)

```

Is this enforces a zero deformation condition only in the y-direction (sub-space 1).

---

<div class="post-metadata">

**Author:** ![Thiago](https://avatars.discourse-cdn.com/v4/letter/t/ec9cab/32.png) [@Thiago](https://fenicsproject.discourse.group/u/Thiago)\
**Post date:** [January 31, 2020, 10:16am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/10 "2020-01-31T10:16:02Z")

</div>

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
*** -------------------------------------------------------------------------
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 31, 2020, 11:53pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/11 "2020-01-31T23:53:20Z")

</div>

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:

```auto
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:

```auto
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

```

---

<div class="post-metadata">

**Author:** ![Thiago](https://avatars.discourse-cdn.com/v4/letter/t/ec9cab/32.png) [@Thiago](https://fenicsproject.discourse.group/u/Thiago)\
**Post date:** [February 3, 2020, 9:56am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/12 "2020-02-03T09:56:55Z")

</div>

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)?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 3, 2020, 10:15am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/13 "2020-02-03T10:15:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![amy](https://avatars.discourse-cdn.com/v4/letter/a/4da419/32.png) [@amy](https://fenicsproject.discourse.group/u/amy)\
**Post date:** [April 8, 2020, 6:37am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/14 "2020-04-08T06:37:16Z")

</div>

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])
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 8, 2020, 7:35am UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/15 "2020-04-08T07:35:35Z")

</div>

You should apply bc\_p to the variational problem, as you do with bc\_base.

---

<div class="post-metadata">

**Author:** ![amy](https://avatars.discourse-cdn.com/v4/letter/a/4da419/32.png) [@amy](https://fenicsproject.discourse.group/u/amy)\
**Post date:** [April 8, 2020, 3:41pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/16 "2020-04-08T15:41:02Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 8, 2020, 4:11pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/17 "2020-04-08T16:11:47Z")

</div>

This is because you were mixing in tolerance inside your subdomain definition, using `near` and `1-tol`. The following code works.

```auto
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])

```

---

<div class="post-metadata">

**Author:** ![amy](https://avatars.discourse-cdn.com/v4/letter/a/4da419/32.png) [@amy](https://fenicsproject.discourse.group/u/amy)\
**Post date:** [April 8, 2020, 4:58pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/18 "2020-04-08T16:58:45Z")

</div>

If I want to have a point deformation rather than a point force in the end of the cantilever, what should I change?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 8, 2020, 5:10pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/19 "2020-04-08T17:10:16Z")

</div>

Right now you have a point deformation. Point force is Added with dirac delta functions.

---

<div class="post-metadata">

**Author:** ![amy](https://avatars.discourse-cdn.com/v4/letter/a/4da419/32.png) [@amy](https://fenicsproject.discourse.group/u/amy)\
**Post date:** [April 8, 2020, 5:25pm UTC](https://fenicsproject.discourse.group/t/help-a-starter/2324/20 "2020-04-08T17:25:48Z")

</div>

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?

[Next page](https://fenicsproject.discourse.group/t/help-a-starter/2324.md?page=2)