# Ferrite 2D problem unphysical results

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.

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

dΩ = getdetJdV(cellvalues, q_point)

for i in 1:getnbasefunctions(cellvalues)
∇ε = tovoigt(SVector,
offdiagscale=2.0)

for j in 1:getnbasefunctions(cellvalues)
∇δε = tovoigt(SVector,
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)
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
cellvalues = CellVectorValues(qr, Lagrange{2,RefCube,1}(), Lagrange{2,RefCube,1}())
facevalues = FaceVectorValues(qr_face, Lagrange{2,RefCube,1}(), Lagrange{2,RefCube,1}())

# DOF and ConstraintHandler
dh = DofHandler(grid)
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()

``````

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

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