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.