BoundsError for J2Plasticity

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 = 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ν)

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

# 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

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

# 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

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)

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

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.

Have you had a look at

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.

1 Like

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

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, so I would suggest you stick to one place in order to not duplicate efforts.)

2 Likes

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

I guess you are referring to BoundsError for Plasticity of material · Discussion #1816 · JuliaLang/www.julialang.org · GitHub, 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.

2 Likes