# BoundsError for J2Plasticity

**URL:** <https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [January 23, 2023, 1:56pm UTC](https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408 "2023-01-23T13:56:29Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Prashanth\_Sudeendran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/prashanth_sudeendran/32/46132_2.png) [@Prashanth\_Sudeendran](https://discourse.julialang.org/u/Prashanth_Sudeendran)\
**Post date:** [January 23, 2023, 1:56pm UTC](https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408/1 "2023-01-23T13:56:30Z")

</div>

Hello all,

I was creating a small variation for a project that I am working for J2Plasticity.  
I ran into this error and I was unable to overcome this and find some possible ways to get out of this situation.  
I have posted screenshots of error and also where the error statements have come from. It would really helpful if someone could find a possible solution for this case.

error being

Info : Meshing 1D…  
Info : [0%] Meshing curve 1 (Line)  
Info : [30%] Meshing curve 2 (Line)  
Info : [50%] Meshing curve 3 (Line)  
Info : [80%] Meshing curve 4 (Line)  
Info : Done meshing 1D (Wall 0s, CPU 0s)  
Info : Meshing 2D…  
Info : Meshing surface 1 (Plane, Frontal-Delaunay)  
Info : Done meshing 2D (Wall 0.00104094s, CPU 0s)  
Info : 76 nodes 154 elements  
Info : Writing ‘C:\Users\prash\AppData\Local\Temp\jl\_GuAloP\mesh.msh’…  
Info : Done writing ‘C:\Users\prash\AppData\Local\Temp\jl\_GuAloP\mesh.msh’  
Info : Reading ‘C:\Users\prash\AppData\Local\Temp\jl\_GuAloP\mesh.msh’…  
Info : 9 entities  
Info : 76 nodes  
Info : 150 elements  
Info : Done reading ‘C:\Users\prash\AppData\Local\Temp\jl\_GuAloP\mesh.msh’

Starting Netwon iterations:

Time step [@time](https://github.com/time) = 1:  
BoundsError: attempt to access 0×0 SparseMatrixCSC{Bool, Int64} with 0 stored entries at index [1, 1]

Stacktrace:  
[1] throw\_boundserror(A::SparseMatrixCSC{Bool, Int64}, I::Tuple{Int64, Int64})  
@ Base .\abstractarray.jl:703  
[2] checkbounds  
@ .\abstractarray.jl:668 [inlined]  
[3] getindex  
@ C:\Users\prash\AppData\Local\Programs\Julia-1.8.2\share\julia\stdlib\v1.8\SparseArrays\src\sparsematrix.jl:2222 [inlined]  
[4] onboundary  
@ C:\Users\prash.julia\packages\Ferrite\Ck78N\src\iterators.jl:90 [inlined]  
[5] assemble\_cell!(Ke::Matrix{Float64}, re::Vector{Float64}, cell::CellIterator{2, Triangle, Float64, DofHandler{2, Float64, Grid{2, Triangle, Float64}}}, cellvalues::CellVectorValues{2, Float64, RefTetrahedron, 4}, facevalues::FaceVectorValues{2, Float64, RefTetrahedron, 4}, grid::Grid{2, Triangle, Float64}, material::J2Plasticity{Float64, SymmetricTensor{4, 2, Float64, 9}}, ue::Vector{Float64}, state::SubArray{MaterialState{Float64, SymmetricTensor{2, 2, Float64, 3}}, 1, Matrix{MaterialState{Float64, SymmetricTensor{2, 2, Float64, 3}}}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, state\_old::SubArray{MaterialState{Float64, SymmetricTensor{2, 2, Float64, 3}}, 1, Matrix{MaterialState{Float64, SymmetricTensor{2, 2, Float64, 3}}}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}, t::Vec{2, Float64})  
@ Main .\In[1]:211  
[6] doassemble(cellvalues::CellVectorValues{2, Float64, RefTetrahedron, 4}, facevalues::FaceVectorValues{2, Float64, RefTetrahedron, 4}, K::SparseMatrixCSC{Float64, Int64}, grid::Grid{2, Triangle, Float64}, dh::DofHandler{2, Float64, Grid{2, Triangle, Float64}}, material::J2Plasticity{Float64, SymmetricTensor{4, 2, Float64, 9}}, u::Vector{Float64}, states::Matrix{MaterialState{Float64, SymmetricTensor{2, 2, Float64, 3}}}, states\_old::Matrix{MaterialState{Float64, SymmetricTensor{2, 2, Float64, 3}}}, t::Vec{2, Float64})  
@ Main .\In[1]:180  
[7] solve()  
@ Main .\In[1]:282  
[8] top-level scope  
@ In[1]:325  
[9] eval  
@ .\boot.jl:368 [inlined]  
[10] include\_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)  
@ Base .\loading.jl:1428

program statements being

using Ferrite, Tensors, SparseArrays, LinearAlgebra, FerriteGmsh, Gmsh, Printf

struct J2Plasticity{T, S \<: SymmetricTensor{4, 2, T}}  
G::T # Shear modulus  
K::T # Bulk modulus  
σ₀::T # Initial yield limit  
H::T # Hardening modulus  
Dᵉ::S # Elastic stiffness tensor  
end;

function J2Plasticity(E, ν, σ₀, H)  
δ(i,j) = i == j ? 1.0 : 0.0 # helper function  
G = E / 2(1 + ν)  
K = E / 3(1 - 2ν)

```julia-auto
Isymdev(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
temp(i,j,k,l) = 2.0G *( 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) + ν/(1.0-2.0ν)*δ(i,j)*δ(k,l))
Dᵉ = SymmetricTensor{4, 2}(temp)
return J2Plasticity(G, K, σ₀, H, Dᵉ)

```

end;

struct MaterialState{T, S \<: SecondOrderTensor{2, T}}

# Store “converged” values

ϵᵖ::S # plastic strain  
σ::S # stress  
k::T # hardening variable  
end

function MaterialState()  
return MaterialState(  
zero(SymmetricTensor{2, 2}),  
zero(SymmetricTensor{2, 2}),  
0.0)  
end

function vonMises(σ)  
s = dev(σ)  
return sqrt(3.0/2.0 \* s ⊡ s)  
end;

function setup\_grid(h=0.05)

# Initialize gmsh

Gmsh.initialize()

```julia-auto
# Add the points
p1 = gmsh.model.geo.add_point(0.0, 0.0, h)
p2 = gmsh.model.geo.add_point(4.0, 0.0, h)
p3 = gmsh.model.geo.add_point(4.0, 3.0, h)
p4 = gmsh.model.geo.add_point(0.0, 3.0, h)

# Add the lines
l1 = gmsh.model.geo.add_line(p1, p2)
l2 = gmsh.model.geo.add_line(p2, p3)
l3 = gmsh.model.geo.add_line(p3, p4)
l4 = gmsh.model.geo.add_line(p4, p1)

# Create the closed curve loop and the surface
loop = gmsh.model.geo.add_curve_loop([l1, l2, l3, l4])
surf = gmsh.model.geo.add_plane_surface([loop])

gmsh.model.geo.synchronize()

# Generate a 2D mesh
gmsh.model.mesh.generate(2)

# Create the physical domains
gmsh.model.add_physical_group(1, [l1], -1, "Γ1")
gmsh.model.add_physical_group(1, [l2], -1, "Γ2")
gmsh.model.add_physical_group(1, [l3], -1, "Γ3")
gmsh.model.add_physical_group(1, [l4], -1, "Γ4")
gmsh.model.add_physical_group(2, [surf])

# Save the mesh, and read back in as a Ferrite Grid
grid = mktempdir() do dir
    path = joinpath(dir, "mesh.msh")
    gmsh.write(path)
    togrid(path)
end

# Finalize the Gmsh library
Gmsh.finalize()

return grid

```

end

function compute\_stress\_tangent(ϵ::SymmetricTensor{2, 2}, material::J2Plasticity, state::MaterialState)

# unpack some material parameters

G = material.G  
H = material.H

```julia-auto
# We use (•)ᵗ to denote *trial*-values
σᵗ = material.Dᵉ ⊡ (ϵ - state.ϵᵖ) # trial-stress
sᵗ = dev(σᵗ) # deviatoric part of trial-stress
J₂ = 0.5 * sᵗ ⊡ sᵗ # second invariant of sᵗ
σᵗₑ = sqrt(3.0*J₂) # effective trial-stress (von Mises stress)
σʸ = material.σ₀ + H * state.k # Previous yield limit

φᵗ = σᵗₑ - σʸ # Trial-value of the yield surface

if φᵗ < 0.0 # elastic loading
    return σᵗ, material.Dᵉ, MaterialState(state.ϵᵖ, σᵗ, state.k)
else # plastic loading
    h = H + 3G
    μ = φᵗ / h # plastic multiplier

    c1 = 1 - 3G * μ / σᵗₑ
    s = c1 * sᵗ # updated deviatoric stress
    σ = s + vol(σᵗ) # updated stress

    # Compute algorithmic tangent stiffness ``D = \frac{\Delta \sigma }{\Delta \epsilon}``
    κ = H * (state.k + μ) # drag stress
    σₑ = material.σ₀ + κ # updated yield surface

    δ(i,j) = i == j ? 1.0 : 0.0
    Isymdev(i,j,k,l) = 0.5*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k)) - 1.0/3.0*δ(i,j)*δ(k,l)
    Q(i,j,k,l) = Isymdev(i,j,k,l) - 3.0 / (2.0*σₑ^2) * s[i,j]*s[k,l]
    b = (3G*μ/σₑ) / (1.0 + 3G*μ/σₑ)

    Dtemp(i,j,k,l) = -2G*b * Q(i,j,k,l) - 9G^2 / (h*σₑ^2) * s[i,j]*s[k,l]
    D = material.Dᵉ + SymmetricTensor{4, 3}(Dtemp)

    # Return new state
    Δϵᵖ = 3/2 * μ / σₑ * s # plastic strain
    ϵᵖ = state.ϵᵖ + Δϵᵖ # plastic strain
    k = state.k + μ # hardening variable
    return σ, D, MaterialState(ϵᵖ, σ, k)
end

```

end

function create\_values(interpolation)

# setup quadrature rules

qr = QuadratureRule{2,RefTetrahedron}(2)  
face\_qr = QuadratureRule{1,RefTetrahedron}(3)

```julia-auto
# create geometric interpolation (use the same as for u)
interpolation_geom = Lagrange{2,RefTetrahedron,1}()

# cell and facevalues for u
cellvalues_u = CellVectorValues(qr, interpolation, interpolation_geom)
facevalues_u = FaceVectorValues(face_qr, interpolation, interpolation_geom)

return cellvalues_u, facevalues_u
return qr, face_qr

```

end;

function create\_dofhandler(grid, interpolation)  
dh = DofHandler(grid)  
dim = 2  
push!(dh, :u, dim, interpolation) # add a displacement field with 2components  
close!(dh)  
return dh  
end

function create\_bc(dh, grid)  
dbcs = ConstraintHandler(dh)

# Clamped on the left side

dofs = [1, 2]  
dbc = Dirichlet(:u, getfaceset(grid, “Γ1”), (x,t) → [0.0, 0.0], dofs)  
add!(dbcs, dbc)  
close!(dbcs)  
return dbcs  
end;

function doassemble(cellvalues::CellVectorValues{dim},  
facevalues::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,  
dh::DofHandler, material::J2Plasticity, u, states, states\_old, t) where {dim}  
r = zeros(ndofs(dh))  
assembler = start\_assemble(K, r)  
nu = getnbasefunctions(cellvalues)  
re = zeros(nu) # element residual vector  
ke = zeros(nu, nu) # element tangent matrix

```julia-auto
for (i, cell) in enumerate(CellIterator(dh))
    fill!(ke, 0)
    fill!(re, 0)
    eldofs = celldofs(cell)
    ue = u[eldofs]
    state = @view states[:, i]
    state_old = @view states_old[:, i]
    assemble_cell!(ke, re, cell, cellvalues, facevalues, grid, material,
                   ue, state, state_old, t)
    assemble!(assembler, eldofs, re, ke)
end
return K, r

```

end

function assemble\_cell!(Ke, re, cell, cellvalues, facevalues, grid, material,  
ue, state, state\_old, t)  
n\_basefuncs = getnbasefunctions(cellvalues)  
reinit!(cellvalues, cell)

```julia-auto
for q_point in 1:getnquadpoints(cellvalues)
    # For each integration point, compute stress and material stiffness
    ϵ = function_symmetric_gradient(cellvalues, q_point, ue) # Total strain
    σ, D, state[q_point] = compute_stress_tangent(ϵ, material, state_old[q_point])

    dΩ = getdetJdV(cellvalues, q_point)
    for i in 1:n_basefuncs
        δϵ = shape_symmetric_gradient(cellvalues, q_point, i)
        re[i] += (δϵ ⊡ σ) * dΩ # add internal force to residual
        for j in 1:i # loop only over lower half
            Δϵ = shape_symmetric_gradient(cellvalues, q_point, j)
            Ke[i, j] += δϵ ⊡ D ⊡ Δϵ * dΩ
        end
    end
end
symmetrize_lower!(Ke)

# Add traction as a negative contribution to the element residual `re`:
for face in 1:nfaces(cell)
    if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "Γ3")
        reinit!(facevalues, cell, face)
        for q_point in 1:getnquadpoints(facevalues)
            dΓ = getdetJdV(facevalues, q_point)
            for i in 1:n_basefuncs
                δu = shape_value(facevalues, q_point, i)
                re[i] -= (δu ⋅ t) * dΓ
            end
        end
    end
end

```

end

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;

function solve()

# Define material parameters

E = 200.0e9 # [Pa]  
H = E/20 # [Pa]  
ν = 0.3 # [-]  
σ₀ = 200e6 # [Pa]  
material = J2Plasticity(E, ν, σ₀, H)  
n\_timesteps = 10  
u\_max = zeros(n\_timesteps)  
traction\_magnitude = 1.e7 \* range(0.5, 1.0, length=n\_timesteps)  
h = 0.05  
grid = setup\_grid(h)  
interpolation = Lagrange{2, RefTetrahedron, 1}() # Linear tet with 2 unknowns/nod  
dh = create\_dofhandler(grid, interpolation) # JuaFEM helper function  
dbcs = create\_bc(dh, grid) # create Dirichlet boundary-conditions

```julia-auto
cellvalues, facevalues = create_values(interpolation)

# Pre-allocate solution vectors, etc.
n_dofs = ndofs(dh) # total number of dofs
u = zeros(n_dofs) # solution vector
Δu = zeros(n_dofs) # displacement correction
r = zeros(n_dofs) # residual
K = create_sparsity_pattern(dh); # tangent stiffness matrix

# Create material states. One array for each cell, where each element is an array of material-
# states - one for each integration point
nqp = getnquadpoints(cellvalues)
states = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)]
states_old = [MaterialState() for _ in 1:nqp, _ in 1:getncells(grid)]

# Newton-Raphson loop
NEWTON_TOL = 1 # 1 N
print("\n Starting Netwon iterations:\n")

for timestep in 1:n_timesteps
    t = timestep # actual time (used for evaluating d-bndc)
    traction = Vec((0.0, traction_magnitude[timestep]))
    newton_itr = -1
    print("\n Time step @time = $timestep:\n")
    update!(dbcs, t) # evaluates the D-bndc at time t
    apply!(u, dbcs) # set the prescribed values in the solution vector

    while true; newton_itr += 1

        if newton_itr > 8
            error("Reached maximum Newton iterations, aborting")
            break
        end
        K, r = doassemble(cellvalues, facevalues, K, grid, dh, material, u,
                         states, states_old, traction);
        norm_r = norm(r[Ferrite.free_dofs(dbcs)])

        print("Iteration: $newton_itr \tresidual: $(@sprintf("%.8f", norm_r))\n")
        if norm_r < NEWTON_TOL
            break
        end

        apply_zero!(K, r, dbcs)
        Δu = Symmetric(K) \ r
        u -= Δu
    end

    # Update the old states with the converged values for next timestep
    states_old .= states

    u_max[timestep] = maximum(abs.(u)) # maximum displacement in current timestep
end

# ## Postprocessing
# Only a vtu-file corrsponding to the last time-step is exported.
#
# The following is a quick (and dirty) way of extracting average cell data for export.
mises_values = zeros(getncells(grid))
κ_values = zeros(getncells(grid))
for (el, cell_states) in enumerate(eachcol(states))
    for state in cell_states
        mises_values[el] += vonMises(state.σ)
        κ_values[el] += state.k*material.H
    end
    mises_values[el] /= length(cell_states) # average von Mises stress
    κ_values[el] /= length(cell_states) # average drag stress
end
vtk_grid("plasticity", dh) do vtkfile
    vtk_point_data(vtkfile, dh, u) # displacement field
    vtk_cell_data(vtkfile, mises_values, "von Mises [Pa]")
    vtk_cell_data(vtkfile, κ_values, "Drag stress [Pa]")
end

return u_max, traction_magnitude

```

end

u\_max, traction\_magnitude = solve();

thanks in advance forum.

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [January 23, 2023, 2:11pm UTC](https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408/2 "2023-01-23T14:11:38Z")

</div>

Have you had a look at

> [@Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757):
>
> Welcome to the Julia Discourse! We are enthusiastic about helping Julia programmers, both beginner and experienced. This public service announcement (PSA) outlines best practices when asking for help. Following these points makes it easier for us to help you and more likely you’ll get a prompt, useful answer. Keywords are highlighted to make it easier to refer to specific points. Choose a descriptive title that captures the key part of your question, eg “plots with multiple axes” instead of …

specifically the advice about an MWE? Your code is \>300 lines, which is unlikely to be minimal and means you’ll find few people if anyone willing to debug this. Try and reduce your example to a few lines of code with dummy data which reproduces the error.

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [January 23, 2023, 4:28pm UTC](https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408/3 "2023-01-23T16:28:11Z")

</div>

Definitely agree with nilshg that you need to minimize your examples in order to get helpful answers. However, I happen to know what the problem is, and it is

> [@Prashanth\_Sudeendran](#):
>
> `onboundary(cell, face)`

which is problematic for externally generated grids. It is just an optimization and you can remove that completely.

(You also asked, and got answers, [here](https://github.com/Ferrite-FEM/Ferrite.jl/discussions/589), so I would suggest you stick to one place in order to not duplicate efforts.)

---

<div class="post-metadata">

**Author:** ![Prashanth\_Sudeendran](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/prashanth_sudeendran/32/46132_2.png) [@Prashanth\_Sudeendran](https://discourse.julialang.org/u/Prashanth_Sudeendran)\
**Post date:** [January 23, 2023, 5:04pm UTC](https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408/4 "2023-01-23T17:04:08Z")

</div>

Thanks for the suggestion , however, I was directed to this page by one other user. So, I asked on this forum.  
I will execute and see how this works. Thanks once again

---

<div class="post-metadata">

**Author:** ![fredrikekre](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikekre/32/1688_2.png) [@fredrikekre](https://discourse.julialang.org/u/fredrikekre)\
**Post date:** [January 23, 2023, 5:07pm UTC](https://discourse.julialang.org/t/boundserror-for-j2plasticity/93408/5 "2023-01-23T17:07:48Z")

</div>

> [@Prashanth\_Sudeendran](#):
>
> Thanks for the suggestion , however, I was directed to this page by one other user. So, I asked on this forum.

I guess you are referring to [BoundsError for Plasticity of material · JuliaLang/www.julialang.org · Discussion #1816 · GitHub](https://github.com/JuliaLang/www.julialang.org/discussions/1816), but that repository is for the Julia webpage and has nothing to do with what you are asking about, which is why you was directed here.

It is perfectly fine to use this forum of course, that is not what I meant. I was just mentioning that you had already gotten some help in the other discussion, and it might make sense to focus on one place, whether here or there.