# Ferrite 2D problem unphysical results
**URL:** https://discourse.julialang.org/t/ferrite-2d-problem-unphysical-results/107037
**Category:** General Usage
**Tags:** pde, numerics, fem
**Created:** [December 2, 2023, 7:46am UTC](https://discourse.julialang.org/t/ferrite-2d-problem-unphysical-results/107037 "2023-12-02T07:46:19Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@c\_sell](https://discourse.julialang.org/u/c_sell)
#### Post date: [December 2, 2023, 7:46am UTC](https://discourse.julialang.org/t/ferrite-2d-problem-unphysical-results/107037/1 "2023-12-02T07:46:19Z")
Hello,
I am using the latest version of Ferrite (Ferrite v0.3.14). I’ve simplified my problem to a small 2D rectangular structure (height = 1; width = 2). The left side is clamped, and on the right-hand side, there is stress in the normal direction. Therefore, I would expect a symmetrical elongation, and not the result shown in the picture below.
Here is the running Ferrite code. I am thankful for any ideas.
```julia
using Ferrite, StaticArrays, SparseArrays, LinearAlgebra
using Plots
using Infiltrator
function compute_u(cellvalues, facevalues,dh,K,f,C,dbc)
Ke = zeros(Float64,8,8)
fe = zeros(Float64,8)
assembler = start_assemble(K,f)
for cell in CellIterator(dh)
Ke .= 0.0
fe .= 0.0
reinit!(cellvalues, cell)
for q_point in getnquadpoints(cellvalues)
dΩ = getdetJdV(cellvalues, q_point)
for i in 1:getnbasefunctions(cellvalues)
∇ε = tovoigt(SVector,
shape_symmetric_gradient(cellvalues, q_point, i),
offdiagscale=2.0)
for j in 1:getnbasefunctions(cellvalues)
∇δε = tovoigt(SVector,
shape_symmetric_gradient(cellvalues, q_point, j),
offdiagscale=2.0)
Ke[i,j] += ∇ε'*C*∇δε*dΩ
end
end
end
for face_id in 1:nfaces(cell)
if onboundary(cell,face_id) && (cellid(cell),face_id) ∈ getfaceset(dh.grid,"pressure")
reinit!(facevalues, cell,face_id)
@show (cellid(cell),face_id)
for q_point in 1:getnquadpoints(facevalues)
dΓ = getdetJdV(facevalues, q_point)
for i in 1:getnbasefunctions(facevalues)
fe[i] += (facevalues.normals[1] ⋅ shape_value(facevalues, q_point, i)) * dΓ
end
end
end
end
assemble!(assembler,celldofs(cell),Ke,fe)
end
apply!(K,f,dbc)
u = Symmetric(K)\f
return u
end
function FEM_MMWE()
## Create Domain ############
corners = [Vec{2}((0.0, 0.0)),
Vec{2}((2.0, 0.0)),
Vec{2}((2.0, 1.0)),
Vec{2}((0.0, 1.0))]
n = 2
grid = generate_grid(Quadrilateral, (n, n), corners);
# node-/facesets for boundary conditions
addnodeset!(grid, "clamped", x -> x[1] ≈ 0.0)
addfaceset!(grid, "pressure", x -> x[1] ≈ 2.0);
## Define Interplotion and QuadratureRule for cells and domains
qr = QuadratureRule{2,RefCube}(3)
cellvalues = CellVectorValues(qr, Lagrange{2,RefCube,1}(), Lagrange{2,RefCube,1}())
qr_face = QuadratureRule{1,RefCube}(3)
facevalues = FaceVectorValues(qr_face, Lagrange{2,RefCube,1}(), Lagrange{2,RefCube,1}())
# DOF and ConstraintHandler
dh = DofHandler(grid)
add!(dh, :u, 2) # displacement
close!(dh)
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getnodeset(dh.grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2]))
close!(dbc)
## Define Material Properties
E = 210e03; ν = 0.3
λ = E/(1-ν^2)
C = Tensor{2,3,Float64}(λ*[1.0 ν 0.0; ν 1.0 0.0; 0.0 0.0 0.5*(1.0-ν)])
K = create_sparsity_pattern(dh)
f = zeros(Float64,size(K,1))
## Compute u
u = compute_u(cellvalues,facevalues,dh,K,f,C,dbc) #compute_u(cellvalues, facevalues,dh,K,f,C,dbc)
### POSTPROCESSING ###
u_nodes = reshape_to_nodes(dh,u,:u)
u_1 = u_nodes[1,:]
u_2 = u_nodes[2,:]
u_coords = getcoordinates.(grid.nodes)
x = getindex.(u_coords,1)
y = getindex.(u_coords,2)
# remove duplicated x coor
x = unique(round.(x,digits = 5))
y = unique(round.(y,digits = 5))
display(x)
display(y)
X = [x_i for x_i in x, _ in y]
Y = [y_i for _ in x, y_i in y]
u_1 = reshape(u_1,length(x),length(y))
u_2 = reshape(u_2,length(x),length(y))
plotly()
surface(X,Y,u_1,title = "x displacement")
display(surface!())
surface(X,Y,u_2,title = "y displacement")
display(surface!())
end
FEM_MMWE()
```
---
### Author:  [@koehlerson](https://discourse.julialang.org/u/koehlerson)
#### Post date: [December 13, 2023, 2:01pm UTC](https://discourse.julialang.org/t/ferrite-2d-problem-unphysical-results/107037/2 "2023-12-13T14:01:08Z")
Hey there,
the problem is that you are unintentionally integrating with one Gauss point the stiffness matrix, due to the missing `1:` in `1:getnquadpoints(cellvalues)`
```diff
17c16
< for q_point in getnquadpoints(cellvalues)
---
> for q_point in 1:getnquadpoints(cellvalues)
```
Further, I’d like to point out that the cast into `tovoigt` can be avoided and you can directly use the entities returned by `shape_symmetric_gradient` if you change the line `Ke[i,j] += ∇ε'*C*∇δε*dΩ` to `Ke[i,j] += ∇ε ⊡ C ⊡ ∇δε*dΩ`. If you want to use it this way, you should change `C` to a `SymmetricTensor{4,2}` and construct it as advertised in the docs of Tensors or Ferrite, see e.g., [Demos · Tensors.jl](https://ferrite-fem.github.io/Tensors.jl/stable/demos/#Creating-the-linear-elasticity-tensor)
---
### Author:  [@termi-official](https://discourse.julialang.org/u/termi-official)
#### Post date: [December 13, 2023, 2:21pm UTC](https://discourse.julialang.org/t/ferrite-2d-problem-unphysical-results/107037/3 "2023-12-13T14:21:47Z")
I would also highly recommend to either use FerriteViz.jl or ParaView for visualization and not trying to do it manually with Plots.jl.
Also note that the Ferrite devs are more active on Slack in #ferrite-fem as well as on Zulip in #Ferrite.jl .