hi all,
im trying to run an hyperplastic problem with 2 different materials and its always fails after few iterations ,
any idea what can i improve?
im not sure which parts of the code to share
this is the defenition of my mesh
E1, E2, nu = 1,10,0.2
mesh = UnitCubeMesh(24, 16, 16)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 1 / 2 + DOLFIN_EPS
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1 / 2 - DOLFIN_EPS
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
x0 = AutoSubDomain(lambda x: near(x[0], 0))
x1 = AutoSubDomain(lambda x: near(x[0], 1))
y0 = AutoSubDomain(lambda x: near(x[1], 0))
z0 = AutoSubDomain(lambda x: near(x[2], 0))
x0.mark(boundary_parts, 1)
y0.mark(boundary_parts, 2)
z0.mark(boundary_parts, 3)
x1.mark(boundary_parts, 4)
class K(fe.UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
E_global = K(materials, E1, E2, degree=0)
DG = FunctionSpace(mesh, "DG", 2)
materials_function = Function(DG)
E_values = [E1, E2]
materials_function = project(E_global, DG)
MU1, LAMBDA1 = Constant(E1/(2*(1 + nu))), Constant(E1*nu/((1 + nu)*(1 - 2*nu)))
MU2, LAMBDA2 = Constant(E2/(2*(1 + nu))), Constant(E2*nu/((1 + nu)*(1 - 2*nu)))
class MUMU(fe.UserExpression):
def __init__(self, materials, MU_0, MU_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.MU_0 = MU_0
self.MU_1 = MU_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.MU_0
else:
values[0] = self.MU_1
def value_shape(self):
return ()
MU = MUMU(materials, MU1, MU2, degree=0)
class lamabada(fe.UserExpression):
def __init__(self, materials, lam_0, lam_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.lam_0 = lam_0
self.lam_1 = lam_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.lam_0
else:
values[0] = self.lam_1
def value_shape(self):
return ()
LAMBDA = lamabada(materials, LAMBDA1, LAMBDA2, degree=0)
of course im solving non linear problem
# Create nonlinear variational problem and solve
problem = NonlinearVariationalProblem(F, _u_p, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.solve()
this is the error i get :
Newton iteration 5: r (abs) = 5.658e+22 (tol = 1.000e-10) r (rel) = 1.210e+20 (tol = 1.000e-06)
Traceback (most recent call last):
File “/home/mirialex/PycharmProjects/hyperelastic_model/new_from_internet.py”, line 228, in
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
*** fenics-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: ubuntu
*** -------------------------------------------------------------------------