hi all,
im facing an issue with my solver with the folllwonig error
Traceback (most recent call last):
File "/home/mirialex/PycharmProjects/fibers/yes_i_can.py", line 425, in <module>
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** https://fenicsproject.discourse.group/
***
*** 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 = meshio.read(f"{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()
dolfin.XDMFFile(f"{mesh_name}_tetra.xdmf").read(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
dolfin_mesh_function.set_values(dolfin_mesh_function_array)
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
dolfin.XDMFFile(f"{mesh_name}.xdmf").write(dolfin_mesh)
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").write(dolfin_subdomains)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").write(dolfin_boundaries)
# 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
os.remove(f"{mesh_name}_tetra.xdmf")
os.remove(f"{mesh_name}_tetra.h5")
os.remove(f"{mesh_name}_triangle.xdmf")
os.remove(f"{mesh_name}_triangle.h5")
mesh_name = "my_perfect_model"
mesh = dolfin.Mesh()
dolfin.XDMFFile(f"{mesh_name}.xdmf").read(mesh)
materials = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
dolfin.XDMFFile(f"{mesh_name}_subdomains.xdmf").read(materials)
print(f"Subdomains: {np.unique(materials.array())}")
boundary_parts = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
dolfin.XDMFFile(f"{mesh_name}_boundaries.xdmf").read(boundary_parts)
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
u.interpolate(initial_guess)
# 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
solver.solve()
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."
try:
# 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(np.dot(defGrad(point), defGrad(point)))][0])
stress_values_all.append(stress123(point)[0])
stress_values_all1.append(np.dot(stress_values_all, weights) * absJ)
deformation_gradient_values_all2.append(np.dot(deformation_gradient_values_all, weights) * absJ)
#stress_values_all=[(np.sqrt(np.dot(stress123(point), stress123(point))))][0]
except Exception as e:
print(f"Could not evaluate at {point}: {e}")
print(stress_values_all1)
print(deformation_gradient_values_all2)
print(stress_values_all)
print(deformation_gradient_values_all)
deformation_gradient_values_all=[]
stress_values_all=[]
# time increment
t += float(dt)
#rest of the code