BoundsError for Plasticity of material

[PrashanthSudeendran]

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.

[logankilpatrick]

Hey! Questions like these should be posts at https://discourse.JuliaLang.org