Solver cannot converge

hi all,

im facing an issue with my solver with the folllwonig error

Traceback (most recent call last):
  File "/home/mirialex/PycharmProjects/fibers/", line 425, in <module>

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset:  ubuntu

i tried to increase the amount of iterations but after 1000 iterations the result starting to give Nan values

i tried to change also my initial guess, to use either (0,0,0) or (x[0],x[1],x[2]) but non of them helped.

i tried to lower the absolute and relative tolerace
i tried to search on previous post but did not found any solution

im posting parts of my code(based on few other codes and help from the community on other posts) that i think relevant , posting only parts so i will be easier to read, let me know if to share other parts as well
also attaching my geo and msh files Dropbox

would appreciate if someone has any advise how to make it to converge

note- it does able to converge if im not using the time loop

import os
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np
import meshio
import dolfin

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

#definnig the young modoulus
E1, E2, nu = 20e6, 20e6, 0.49

mesh_name = "my_perfect_model"

# Then, read it in through meshio
meshio_mesh ="{mesh_name}.msh")

# Extract mesh corresponding to each cell type, and write it out with meshio
def extract_mesh_by_cell_type(meshio_mesh, cell_type):
    cells = meshio_mesh.get_cells_type(cell_type)
    cell_data = meshio_mesh.get_cell_data("gmsh:physical", cell_type)
    return meshio.Mesh(points=meshio_mesh.points, cells={cell_type: cells}, cell_data={"markers": [cell_data]})

tetra_mesh = extract_mesh_by_cell_type(meshio_mesh, "tetra")
meshio.write(f"{mesh_name}_tetra.xdmf", tetra_mesh)
del tetra_mesh
triangle_mesh = extract_mesh_by_cell_type(meshio_mesh, "triangle")
meshio.write(f"{mesh_name}_triangle.xdmf", triangle_mesh)
del triangle_mesh

del meshio_mesh

# Read back in the meshio mesh from XDMF in order to convert to a dolfin mesh
dolfin_mesh = dolfin.Mesh()

# Read back in the markers in order to attach them to the dolfin mesh
def read_back_in_markers(meshio_xdmf_file, dolfin_mesh):
    dolfin_markers = dolfin.MeshValueCollection("size_t", dolfin_mesh)
    dolfin.XDMFFile(meshio_xdmf_file).read(dolfin_markers, "markers")
    dolfin_mesh_function = dolfin.cpp.mesh.MeshFunctionSizet(dolfin_mesh, dolfin_markers)

    dolfin_mesh_function_array = dolfin_mesh_function.array()
    defined_markers = np.array(list(set(dolfin_markers.values().values())), dtype=dolfin_mesh_function_array.dtype)
    defined_entries = np.isin(dolfin_mesh_function_array, defined_markers)
    dolfin_mesh_function_array[~defined_entries] = 0
    return dolfin_mesh_function

dolfin_subdomains = read_back_in_markers(f"{mesh_name}_tetra.xdmf", dolfin_mesh)
dolfin_boundaries = read_back_in_markers(f"{mesh_name}_triangle.xdmf", dolfin_mesh)

# Now export them to file

# Please now open this three files in paraview and check that the labels are as expected

# You may remove *_tetra.* and *_triangle.*, they are not needed anymore

mesh_name = "my_perfect_model"

mesh = dolfin.Mesh()

materials = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
print(f"Subdomains: {np.unique(materials.array())}")

boundary_parts = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
print(f"Boundaries: {np.unique(boundary_parts.array())}")

V = VectorFunctionSpace(mesh, "Lagrange", 1)

bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 19)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 21)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 23)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 24)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('((1/(time_+1))-1)*0.2'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 20)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 22)

bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
subdomain_ids = [20, 22]
ds = Measure("ds", subdomain_data=boundary_parts, subdomain_id=(20,22))
ds = ds(degree=4)

#rest of the code

du = TrialFunction(V)  # Incremental displacement
v = TestFunction(V)  # Test function
u = Function(V)  # Displacement from previous iteration
B = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T = Constant((0.0, 0.0, 0.0))  # Traction force on the boundary

# Project or interpolate the initial guess onto u
initial_guess = Expression(("x[0]", "x[1]", "x[2]"), degree=1)

# Project or interpolate the initial guess onto u

# Kinematics
d = len(u)
I = Identity(d)  # Identity tensor
F = variable(I + grad(u))  # Deformation gradient
C = F.T * F  # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Stored strain energy density (compressible neo-Hookean model)

W = MU / 2 * (tr(C) - 3 - 2 * ln(J)) + LAMBDA / 2 * pow((J - 1), 2)

n = FacetNormal(mesh)
P = diff(W, F)
# Total potential energy
L = inner(P, grad(v)) * dx - dot(B, v) * dx - dot(T, v) * ds

Je = derivative(L, u, du)

problem = NonlinearVariationalProblem(L, u, bcs=bcs,J=Je)

solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-3
solver.parameters['newton_solver']['absolute_tolerance'] = 1e-2
solver.parameters['newton_solver']['linear_solver'] = 'lu'
solver.parameters['newton_solver']['maximum_iterations'] = 1200
solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'
solver.parameters['newton_solver']['relaxation_parameter'] = 0.8

#rest of the code

while t <= T:
    print('time: ', t)

    # Increase traction
    TT.t = t
    tDirBC.time_ = t
    tDirBC4.time_ = t

    # solving the problem with time dependensy
    F = variable(I + grad(u))
    J = det(F)
    CL = F * F.T
    # Deviatoric part
    dev_stress = MU * (CL - I)

    # Volumetric part
    vol_stress = LAMBDA * (J - 1) * I

    sigma = (1 / J) * dev_stress + vol_stress
    u.rename("u", "displacement")

    DF = I + grad(u)  # Compute the deformation gradient
    defGrad.assign(project(DF, V_tensor))  # Project the deformation gradient
    sigma_xx = sigma # Extravgacting the xx component of the stress tensor
    stress_xx_projected = project(sigma_xx, V_tensor)  # Projecting this scalar component
    stress123.assign(stress_xx_projected)  # Assigning the projected stress
    # get stretch at a point for plotting
    #trying to create avarage stress and deformation

    for point in physical_points:
        # is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."


            # is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."

            deformation_gradient_values_all.append(defGrad(point)[0])# Store the value at the current point
            #deformation_gradient_values_all=([np.sqrt(, defGrad(point)))][0])
            stress_values_all1.append(, weights) * absJ)
            deformation_gradient_values_all2.append(, weights) * absJ)

            #stress_values_all=[(np.sqrt(, stress123(point))))][0]

        except Exception as e:
            print(f"Could not evaluate at {point}: {e}")



    # time increment
    t += float(dt)

#rest of the code

It is hard to give any specific guidance without having the mesh etc.
In general, try first to switch to a small strain formulation (just changing the potential energy) to see if results make sense in terms of bcs and loading
Then go back to the finite strain setting and start with small load steps. In general a correct Newton method should either not converge if you are too far away from the solution or converge in less than 10 iterations

i attached the mesh here

hi @bleyerj ,

i was wondering if you had time to check my mesh that i uploaded.
maybe it can help more to understand the issue.
i did tried to lower the load steps but unfortunately it didnt helped.
before i used this mesh i used a regular build in meshbox

length, width, height = 1.0, 0.2, 0.2  # Dimensions of the rod
num_elements = 40  # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 4, 4)

and its worked good, this issue of not converging started when i started using this mesh that i uploaded to the dropbox in the previous comment