# No convergence for outsourced function which is included in the variation formulation

Hello everyone,

i am currently writing my first python package for modeling a hyperelastic problem. Therefor I want to write the different material laws and kinematic quantities in separate files and call them from a central file.
Unfortunately, the solution does not converge if the ufl variables for the kinematic quantities are not defined directly at the weak form of the problem itself.

Since the problem cannot be easily translated into a short MWE, I’ll try the following example.
I assume that I have fundamentally misunderstood or misapplied something here and the error can be easily detected by someone more experienced than me.

Code that works:

``````# FILE A: material_models
import ufl
from dolfinx import fem
from petsc4py.PETSc import ScalarType

def neo_hooke(domain, u):
I = ufl.variable(ufl.Identity(len(u)))
C = ufl.variable(F.T * F)
I_1 = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu / 2) * (I_1 - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
P = ufl.diff(psi, F)
S = 2 * ufl.diff(psi, C)
return P, S
``````

Code that does not work:

``````# FILE A: material_models

import ufl
from dolfinx import fem
from petsc4py.PETSc import ScalarType
from kinematics import *

def neo_hooke(domain, u):
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu / 2) * (I1(u) - 3) - mu * ufl.ln(jacobian_defgrad(u)) + (lmbda / 2) * (ufl.ln(jacobian_defgrad(u)))**2
S = 2 * ufl.diff(psi, right_cauchygreen(u))
return P, S

# FILE B: kinematics
import ufl
import numpy as np

def identity_tensor(u):
d = len(u)
return ufl.variable(ufl.Identity(d))

I = identity_tensor(u)

return ufl.variable(ufl.det(F))

def right_cauchygreen(u):
return ufl.variable(F.T * F)

def I1(u):
C = right_cauchygreen(u)
I1 = ufl.tr(C)
return ufl.variable(I1)
``````

Why cant you make this into an MWE?

1. Simply create a unit square or cube
2. Define the function space `V` and the function `u`.
3. Create the variational form in 2 ways.

This follows here:

``````import ufl
from dolfinx import fem, mesh
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
def neo_hooke(domain, u):
I = ufl.variable(ufl.Identity(len(u)))
C = ufl.variable(F.T * F)
I_1 = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu / 2) * (I_1 - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
P = ufl.diff(psi, F)
S = 2 * ufl.diff(psi, C)
return P, S

def neo_hooke2(domain, u):
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
psi = (mu / 2) * (I1(u) - 3) - mu * ufl.ln(jacobian_defgrad(u)) + (lmbda / 2) * (ufl.ln(jacobian_defgrad(u)))**2
S = 2 * ufl.diff(psi, right_cauchygreen(u))
return P, S

def identity_tensor(u):
d = len(u)
return ufl.variable(ufl.Identity(d))

I = identity_tensor(u)

return ufl.variable(ufl.det(F))

def right_cauchygreen(u):
return ufl.variable(F.T * F)

def I1(u):
C = right_cauchygreen(u)
I1 = ufl.tr(C)
return ufl.variable(I1)

domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = fem.VectorFunctionSpace(domain, ("Lagrange", 1))
u = fem.Function(V)

neo_hooke(domain, u)

neo_hooke2(domain, u)
``````

Neither of these produce an error for me, so there is no way of debugging your issue.

MWE here we go:

FILE 1: test.py

``````import numpy as np
import ufl
from dolfinx import fem, mesh, log
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from dolfinx import nls
from kinematix import *

# domain geometries
length = 0.25
height = 0.02

# dimensions of indenter
D = 0.15
indentation_depth = height/2.0
center_x = length / 2
center_y = (height + (D/2)) - indentation_depth

domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0],[length, height]], [40, 9])
x = ufl.SpatialCoordinate(domain)
n = ufl.FacetNormal(domain)

def bottom(x):
return np.isclose(x, 0)

def top(x):
return np.isclose(x, height)

tdim = domain.topology.dim
fdim = tdim - 1

bottom_facets = mesh.locate_entities_boundary(domain,fdim,bottom)
top_facets = mesh.locate_entities_boundary(domain,fdim,top)

marked_facets = np.hstack([bottom_facets, top_facets])
marked_values = np.hstack([np.full_like(bottom_facets, 1), np.full_like(top_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

V = fem.VectorFunctionSpace(domain, ("CG",1))

u = fem.Function(V, name="Displacement")
du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)

u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_tag.find(1)), V)

# definition of gap function
def gapp(u, x):
R = ufl.sqrt(pow(((x - center_x) + u), 2) + pow(((x - center_y) + u), 2))
return  R - (D/2)

# maculay bracket
def maculay(x):
return (x+abs(x))/2

# kinematic quantities
I = ufl.variable(ufl.Identity(len(u)))
C = ufl.variable(F.T * F)
I_1 = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Neo-Hooke
E = ScalarType(1.0e4)
nu = ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))

################################## questionable part ##################################
psi = (mu / 2) * (I_1 - 3)   - mu * ufl.ln(J)                   + (lmbda / 2) * (ufl.ln(J))**2 # works
#psi = (mu / 2) * (I1(u) - 3) - mu * ufl.ln(jacobian_defgrad(u)) + (lmbda / 2) * (ufl.ln(jacobian_defgrad(u)))**2 # does not work
#######################################################################################

# 1st Piola-Kirchhoff
P = ufl.diff(psi, F)

# Penalty Stiffness
kpen = fem.Constant(domain, ScalarType(1e7))

# penalty term for contact
penalty = kpen * n *maculay(-gapp(x, u))

F2 = kpen * ufl.inner(u_, n*maculay(-gapp(x, u)))*ds(2)

form = F1 + F2
Jac = ufl.derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=Jac)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

solver.atol = 1e-6
solver.rtol = 1e-6
solver.convergence_criterion = "incremental"
solver.report = True
ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)

w, converged = solver.solve(u)
assert(converged)
print(f"Number of interations: {w:d}")

with XDMFFile(domain.comm, "mwe_package/result/output.xdmf", "w") as xdmf:
xdmf.write_mesh(domain)
xdmf.write_function(u)
``````

FILE 2: kinematix.py

``````import ufl

def identity_tensor(u):
d = len(u)
return ufl.variable(ufl.Identity(d))

I = identity_tensor(u)

return ufl.variable(ufl.det(F))

def right_cauchygreen(u):