# Compression of a box

**URL:** <https://fenicsproject.discourse.group/t/compression-of-a-box/6676>\
**Category:** Uncategorized\
**Created:** [September 27, 2021, 10:48am UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676 "2021-09-27T10:48:30Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 27, 2021, 10:48am UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/1 "2021-09-27T10:48:30Z")

</div>

Hi,  
I am a beginner in fenics, i would like to simulate the compression of a box, i did the simulation but resulats (using paraview) was illogical especially for the Von Mises. Can you find the mistakes made in my code:

```auto
from fenics import *
from ufl import nabla_div
rho = 1500e-9
E = 3000.2
nu = 0.38
mu = E/(2*(1+nu))
lambda_ = (nu*E)/((1+nu)*(1-2*nu))
tol=1E-9
mesh = BoxMesh(Point(0,0,0),Point(100,100,100),20,20,20)
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]<0+tol
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]>100-tol
sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryBottom.mark(sub_domains, 1)
boundaryTop = BoundaryTop()
boundaryTop.mark(sub_domains, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V, Constant((0,0,-30)), boundaryTop)
bc=[bc1,bc2]
d = u.geometric_dimension()
def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return 2*mu*epsilon(u) + lambda_*nabla_div(u)*Identity(d)
a = inner(sigma(u), epsilon(v))*dx
f= Constant((0,0,0))
L = dot(f,v)*ds
u=Function(V)
solve(a==L, u, bc)
vtkfile = File('u.pvd')
vtkfile << u
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
W = FunctionSpace(mesh, 'CG', 0)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W)
vtkfile = File('s.pvd')
vtkfile << von_Mises

```

 ![mises](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/993c5a12764c1b75984719e8ce3516dec65937e7.jpeg)

Thanks in advance for your help!

---

<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:** [September 27, 2021, 11:47am UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/2 "2021-09-27T11:47:36Z")

</div>

Please format your code using 3x` such that it is possible to copy-paste the code to an editor and run it without modifications.  
Please also note that your VonMises stresses should be projected into a DG 0 space, as the displacement lives in CG 1. This has also been noted in other questions such as:

> [@Results obtained by FEniCS mismatches Abaqus](https://fenicsproject.discourse.group/t/results-obtained-by-fenics-mismatches-abaqus/3268/2):
>
> The code you supplied is not complete, and therefore there is not much one can do to investigate this. Note that the von Mises stresses should be projected into W (the discontinuous space).

---

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 27, 2021, 12:26pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/3 "2021-09-27T12:26:51Z")

</div>

I think that now the code can be copied, Thanks in advance for your help! But beside if I project the VonMises stresses into a DG 0 space, i have the same distributions and shape

---

<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:** [September 27, 2021, 12:46pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/4 "2021-09-27T12:46:03Z")

</div>

As long as you save it to a PVD, the data will be stored at vertices, and not be properly discontinuous, you should use `XDMFile.write_checkpoint`, as shown in [Loading xdmf data back in - #4 by dokken](https://fenicsproject.discourse.group/t/loading-xdmf-data-back-in/1925/4)

---

<div class="post-metadata">

**Author:** ![bleyerj](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bleyerj/32/6020_2.png) [@bleyerj](https://fenicsproject.discourse.group/u/bleyerj)\
**Post date:** [September 27, 2021, 1:45pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/5 "2021-09-27T13:45:51Z")

</div>

Hi,  
if you restrain the lateral displacements on the top and bottom surfaces, you are not simulating simple compression. This will result in high stress concentrations near the border/corners of you domain as illustrated in your figure.  
Consider something like:

```auto
bc1= DirichletBC(V.sub(2), Constant(0.), boundaryBottom)
bc2= DirichletBC(V.sub(2), Constant(-30), boundaryTop)

```

You also need to remove rigid body motions by fixing horizontal displacement and rotations about the Z axis : either by enforcing u\_x=0 (resp. u\_y=0) on a symmetry plane such as x=0 (resp. y=0) or by fixing appropriate displacements at some specific points if you are not satisfied with symmetry conditions

---

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 27, 2021, 4:21pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/6 "2021-09-27T16:21:37Z")

</div>

I saved the Von mises in XDMF file using:

```auto
save = XDMFFile('u_xdmf.xdmf')
save.write_checkpoint(von_Mises, "mises", 0)

```

But I get the same result as before, the distribution of von mises if different than the resulat that i get using comsol 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:** [September 27, 2021, 4:34pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/7 "2021-09-27T16:34:17Z")

</div>

see @kamensky’s comment above

---

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 28, 2021, 11:18am UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/8 "2021-09-28T11:18:12Z")

</div>

Thank you very much for the help. Now, the von mises distribution is logical.

 ![70](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/2819279777ffc0b88c053388a24b5b1733b7cc7e.png)

but there is a difference between the comsol results and fenics results, so i decided to use the same mesh in both, I exported the mesh created in comsol in stl format then i converted into xml format using meshio to use it in fenics. But, using the new mesh leads to illogical results.

 ![90](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/6/6618134c7c7fe8f726f32ab5933f055025183ef6.png)

Have you an idea about the cause and how can I resolve it?  
Thanks in advance for your help!

---

<div class="post-metadata">

**Author:** ![Ekrem\_Ekici](https://avatars.discourse-cdn.com/v4/letter/e/977dab/32.png) [@Ekrem\_Ekici](https://fenicsproject.discourse.group/u/Ekrem_Ekici)\
**Post date:** [September 28, 2021, 2:23pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/9 "2021-09-28T14:23:28Z")

</div>

It would be nice if you rescale the colormaps for sake of comparison.

---

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 28, 2021, 2:29pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/10 "2021-09-28T14:29:29Z")

</div>

It’s done, I changed the color and data range

---

<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:** [September 28, 2021, 4:12pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/11 "2021-09-28T16:12:39Z")

</div>

Please add the code that you are currently executing as to reflect the changes that you have made since the initial post.

---

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 28, 2021, 4:30pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/12 "2021-09-28T16:30:06Z")

</div>

```auto
from fenics import *
from ufl import nabla_div
E = 2500
nu = 0.4
mu = E/(2*(1+nu))
lambda_ = (nu*E)/((1+nu)*(1-2*nu))
tol=1E-6
mesh = Mesh('Box.xml')
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]<0+tol
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]>100-tol
sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryBottom.mark(sub_domains, 1)
boundaryTop = BoundaryTop()
boundaryTop.mark(sub_domains, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V.sub(2), Constant(-30), boundaryTop)
bc=[bc1,bc2]
d = u.geometric_dimension()
def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return 2*mu*epsilon(u) + lambda_*nabla_div(u)*Identity(d)
a = inner(sigma(u), epsilon(v))*dx
f= Constant((0,0,0))
L = dot(f,v)*ds
u=Function(V)
solve(a==L, u, bc)
vtkfile = File('u.pvd')
vtkfile << u
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
W = FunctionSpace(mesh, 'P', 1)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W)
vtkfile = File('s.pvd')
vtkfile << von_Mises
save = XDMFFile('Von_mises.xdmf') 
save.write_checkpoint(von_Mises, "mises", 0) 

```

I has the same illogical distribution when I I project the VonMises stresses into a DG 0 space.  
In this code, as I said before, i converted the box from STL format using meshio into xml format to use it in fenics. When I used the built-in mesh i has logical results, but it’s not the case with an imported mesh

---

<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:** [September 28, 2021, 4:47pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/13 "2021-09-28T16:47:30Z")

</div>

You get quite different results if you use CG 1 or DG 0 visualization. Here is your code with a built in mesh, using DG-0:

```python
**strong text** from fenics import *
from ufl import nabla_div
E = 2500
nu = 0.4
mu = E/(2*(1+nu))
lambda_ = (nu*E)/((1+nu)*(1-2*nu))
tol=1E-6
N = 20
mesh = BoxMesh(Point(0,0,0),Point(100,100,100), N, N, N)
class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]<0+tol
class BoundaryTop(SubDomain):
    def inside(self, x, on_boundary):
        return x[2]>100-tol
sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim())
sub_domains.set_all(0)
boundaryBottom = BoundaryBottom()
boundaryBottom.mark(sub_domains, 1)
boundaryTop = BoundaryTop()
boundaryTop.mark(sub_domains, 2)
V = VectorFunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
bc1= DirichletBC(V, Constant((0,0,0)), boundaryBottom)
bc2= DirichletBC(V.sub(2), Constant(-30), boundaryTop)
bc=[bc1,bc2]
d = u.geometric_dimension()
def epsilon(u):
    return sym(grad(u))
def sigma(u):
    return 2*mu*epsilon(u) + lambda_*nabla_div(u)*Identity(d)
a = inner(sigma(u), epsilon(v))*dx
f= Constant((0,0,0))
L = dot(f,v)*ds
u=Function(V)
solve(a==L, u, bc)
vtkfile = File('u.pvd')
vtkfile << u
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)
W = FunctionSpace(mesh, 'DG', 0)
von_Mises = project(sqrt (3/2.* inner (s, s) ), W)
save = XDMFFile('Von_mises.xdmf') 
save.write_checkpoint(von_Mises, "mises", 0) 

```

and you get:

 ![mises_dg0](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/3/392ffcdc2bfdb5d840137365a975089a776b8af1.png)  
as compared to  
 ![mises_cg1](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/c/c9bb3bde8db158a9ef4dd40d36d0d81bb74dd2a4.png)  
which I claim is radically different solutions, as the magnitude of the stress is different  
DG 0 ranges from (630.882 to 1477.44) on the surface  
CG 1 ranges from (633.903 to 1403.21) on the surface

> [@Salim](#):
>
> When I used the built-in mesh i has logical results, but it’s not the case with an imported mesh

You need to supply your mesh if you want any further guidance on this subject. Could you please provide your original mesh file?

---

<div class="post-metadata">

**Author:** ![Salim](https://avatars.discourse-cdn.com/v4/letter/s/94ad74/32.png) [@Salim](https://fenicsproject.discourse.group/u/Salim)\
**Post date:** [September 28, 2021, 5:01pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/14 "2021-09-28T17:01:29Z")

</div>

[Mesh (Stl format)](https://drive.google.com/file/d/1E8wGuDInR3qCdiaUsXIybbGs4pbzy7od/view?usp=sharing)  
[Mesh (Xml format after the conversion)](https://drive.google.com/file/d/13zk5GmadKs8cLNkx-3yfhLrB-iazEM5K/view?usp=sharing)

---

<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:** [September 28, 2021, 5:41pm UTC](https://fenicsproject.discourse.group/t/compression-of-a-box/6676/15 "2021-09-28T17:41:44Z")

</div>

As you can see from the xml data in your file, your mesh in only a surface mesh (it does not contain tetrahedron, only triangles). This is because your step file only contains the triangular elements.