Mixed elements can be used to overcome locking when the material becomes incompressible. However, for the elements to be stable, they need to fulfil the LBB condition. We here show what happens with a linear / linear displacement pressure element (which does not fulfil the LBB condition). In the numerical example, we consider the Cook's Membrane problem with an applied traction on the right hand side.
In [ ]:
using JuAFEM
using BlockArrays
In [ ]:
function create_cook_grid(nx, ny)
dim = 2
corners = [Vec{dim}((0.0, 0.0)),
Vec{dim}((48.0, 44.0)),
Vec{dim}((48.0, 60.0)),
Vec{dim}((0.0, 44.0))]
grid = generate_grid(Triangle, (nx, ny), corners);
# Extract the left boundary
addnodeset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
return grid
end;
In [ ]:
dim = 2
# Interpolations
ip_u = Lagrange{dim, RefTetrahedron, 1}()
ip_p = Lagrange{dim, RefTetrahedron, 1}()
# Quadrature rules
qr = QuadratureRule{dim , RefTetrahedron}(3)
qr_face = QuadratureRule{dim-1, RefTetrahedron}(3)
cellvalues_u = CellVectorValues(qr, ip_u);
facevalues_u = FaceVectorValues(qr_face, ip_u);
cellvalues_p = CellScalarValues(qr, ip_p);
In [ ]:
# DofHandler
function create_dofhandler(grid)
dh = DofHandler(grid)
push!(dh, :u, dim) # Add a displacement field
push!(dh, :p, 1) # Add a pressure field
close!(dh)
return dh
end;
In [ ]:
# Boundaryconditions
function create_boundaryconditions(dh, grid)
dbc = DirichletBoundaryConditions(dh)
# Add a homogenoush boundary condition on the "clamped" edge
add!(dbc, :u, getnodeset(grid, "clamped"), (x,t) -> zero(Vec{2}), [1,2])
close!(dbc)
t = 0.0
update!(dbc, t)
return dbc
end;
In [ ]:
function symmetrize_lower!(K)
for i in 1:size(K,1)
for j in i+1:size(K,1)
K[i,j] = K[j,i]
end
end
end;
In [ ]:
immutable LinearElasticity{T}
G::T
K::T
end
In [ ]:
function doassemble{dim}(cellvalues_u::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim},
facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
dh::DofHandler, mp::LinearElasticity)
f = zeros(ndofs(dh))
assembler = start_assemble(K, f)
assembler2 = start_assemble(K, f)
nu = getnbasefunctions(cellvalues_u)
np = getnbasefunctions(cellvalues_p)
global_dofs = zeros(Int, nu + np)
fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # Local force vector
Ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # Local stiffness mastrix
t = Vec{2}((0.0, 1/16)) # Traction vector
ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:getnbasefunctions(cellvalues_u)]
for cell in CellIterator(dh)
fill!(Ke, 0)
fill!(fe, 0)
assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)
celldofs!(global_dofs, cell)
assemble!(assembler, global_dofs, fe, Ke)
end
return K, f
end;
In [ ]:
function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)
n_basefuncs_u = getnbasefunctions(cellvalues_u)
n_basefuncs_p = getnbasefunctions(cellvalues_p)
u▄, p▄ = 1, 2
reinit!(cellvalues_u, cell)
reinit!(cellvalues_p, cell)
# We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
@inbounds for q_point in 1:getnquadpoints(cellvalues_u)
for i in 1:n_basefuncs_u
ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))
end
dΩ = getdetJdV(cellvalues_u, q_point)
for i in 1:n_basefuncs_u
divδu = shape_divergence(cellvalues_u, q_point, i)
δu = shape_value(cellvalues_u, q_point, i)
for j in 1:i
Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ
end
end
for i in 1:n_basefuncs_p
δp = shape_value(cellvalues_p, q_point, i)
for j in 1:n_basefuncs_u
divδu = shape_divergence(cellvalues_u, q_point, j)
Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
end
for j in 1:i
p = shape_value(cellvalues_p, q_point, j)
Ke[BlockIndex((p▄, p▄), (i, j))] += - 1/mp.K * δp * p * dΩ
end
end
end
symmetrize_lower!(Ke)
@inbounds for face in 1:nfaces(cell)
if onboundary(cell, face) && (JuAFEM.cellid(cell), face) ∈ getfaceset(grid, "right")
reinit!(facevalues_u, cell, face)
for q_point in 1:getnquadpoints(facevalues_u)
dΓ = getdetJdV(facevalues_u, q_point)
for i in 1:n_basefuncs_u
δu = shape_value(facevalues_u, q_point, i)
fe[i] += (δu ⋅ t) * dΓ
end
end
end
end
end;
In [ ]:
function solve(ν, doexport = true)
# Material
Emod = 1.
Gmod = Emod / 2(1 + ν)
Kmod = Emod * ν / ((1+ν) * (1-2ν))
mp = LinearElasticity(Gmod, Kmod)
# Grid, dofhandler, boundary condition
n = 50
grid = create_cook_grid(n, n)
dh = create_dofhandler(grid)
dbc = create_boundaryconditions(dh, grid)
# Assembly and solve
K = create_sparsity_pattern(dh);
@time K, f = doassemble(cellvalues_u, cellvalues_p, facevalues_u, K, grid, dh, mp);
apply!(K, f, dbc)
u = Symmetric(K) \ f;
# Export
if doexport
vtkfile = vtk_grid("up_$ν", dh, u)
vtk_save(vtkfile)
end
return u
end
for ν in [0.3, 0.4999999]
solve(ν)
end
In [ ]:
u = solve(0.3, false)
Base.Test.@test maximum(u) ≈ 26.13381519901358
println("Cook passed!")
In [ ]:
#= TODO: Mini element
immutable MiniDisplacements{dim, shape, order} <: Interpolation{dim, shape, order} end
JuAFEM.getnbasefunctions(::MiniDisplacements{2, RefTetrahedron, 1}) = 4
function JuAFEM.value!(ip::MiniDisplacements{2, RefTetrahedron, 1}, N::AbstractVector, ξ::Vec{2})
@assert length(N) == 4
JuAFEM.value!(Lagrange{2, RefTetrahedron, 1}(), view(N, 1:3), ξ)
N[4] = N[1] * N[2] * N[3]
return N
end
function JuAFEM.derivative!{T}(ip::MiniDisplacements{2, RefTetrahedron, 1}, dN::AbstractVector{Vec{2, T}}, ξ::Vec{2, T})
@assert length(dN) == 4
ξx, ξy = ξ[1], ξ[2]
JuAFEM.derivative!(Lagrange{2, RefTetrahedron, 1}(), view(dN, 1:3), ξ)
dN[4] = Vec{2, T}((ξy * (1 - 2ξx - ξy),
ξx * (1 - 2ξy - ξx)))
return dN
end
cellvalues_u_mini = CellVectorValues(qr_mini, ip_u_mini);
facevalues_u_mini = FaceVectorValues(qr_face_mini, ip_u_mini);
ip_u_mini = MiniDisplacements{dim, RefTetrahedron, 1}()
qr_mini = QuadratureRule{dim , RefTetrahedron}(3)
qr_face_mini = QuadratureRule{dim-1, RefTetrahedron}(3)
# Integrates along the right boundary
function integrate_gamma(u, facevalues_u, grid, dh)
global_dofs = zeros(Int, ndofs_per_cell(dh))
u_integrated = zero(Vec{2})
for cell in CellIterator(dh)
celldofs!(global_dofs, cell)
up_nodes = u[global_dofs]
u_nodes = up_nodes[1:getnbasefunctions(facevalues_u)]
for face in 1:nfaces(cell)
if onboundary(cell, face) && (JuAFEM.cellid(cell), face) ∈ getfaceset(grid, "right")
reinit!(facevalues_u, cell, face)
for q_point in 1:getnquadpoints(facevalues_u)
dΓ = getdetJdV(facevalues_u, q_point)
u_cell = function_value(facevalues_u, q_point, u_nodes)
u_integrated += u_cell * dΓ
end
end
end
end
return u_integrated
end
=#