Problems with ufl_replace for NonlinearProblem

dolfinx
1

Hi

Working with a non-standard nonlinear solid mechanics problem (so-called change of director field as additional unknown besides displacement field), I am trying to obtain the forms of two stress tensors in the reference/undeformed configuration (P0mf and d0mf) using ufl.replace(), so that I can subsequently subtract these forms from the actual stress tensors forms to ensure that they are both zero tensors in the reference/undeformed configuration (Pmf = Pmf - P0mf and dmf = dmf - d0mf). However, the subtracting seems to remove the stress tensors and their contributions in the variational principle altogether when initializing the NonlinearProblem and solving it, the change of director field unknown is zero regardless. I have implemented the exact same in my C++ code where there is no issue, i.e. the change of director field is not zero as expected.

########################################################################################################################
# Function space
P1 = ufl.VectorElement("CG", msh.ufl_cell(), 2, dim=3)  # Lagrange family, 2nd order
P2 = ufl.VectorElement("CG", msh.ufl_cell(), 2, dim=3)  # Lagrange family, 1st order
VY = fem.FunctionSpace(msh, ufl.MixedElement([P1, P2]))

v, yf = ufl.TestFunctions(VY)     # test functions
uwf = fem.Function(VY)            # current displacement and change of director
u, wf = ufl.split(uwf)

########################################################################################################################
# Kinematics

d = len(u)

Af = fem.Constant(msh, PETSc.ScalarType((1.0/ufl.sqrt(2), 1.0/ufl.sqrt(2), 0.0)))

def director(wf_h):
    """Return an expression for the director field"""
    return Af + wf_h

def defGradTens(u_h):
    """Return an expression for the deformation gradient field"""
    return ufl.Identity(d) + ufl.grad(u_h)

F = ufl.variable(defGradTens(u)) # deformation gradient must be variable to derive 1st Piola-Kirchhoff stress tensor
f = ufl.inv(F) # inverse of deformation gradient
C = F.T * F # right Cauchy-Green deformation tensor
c = f.T * f # Cauchy deformation tensor

# Matrix strain invariants
I1 = ufl.tr(C)
J  = ufl.det(F) # Jacobian

# Fibre strain invariants
af = ufl.variable(director(wf)) # director must be variable to derive director stress vector
L4f = ufl.inner(af, af)

# Fibre-matrix interaction strain invariants
J4f = ufl.inner(c, ufl.outer(af, af))


########################################################################################################################
# Constitutive law

# Set the elasticity parameters

# matrix
E = PETSc.ScalarType(1.0e0)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(msh, E/(2*(1 + nu))) 
lmbda = fem.Constant(msh, E*nu/((1 + nu)*(1 - 2*nu)))

# matrix-fibre bond stiffness
b4Ff = fem.Constant(msh, 1.0e1)  

# fibre stiffness
c4f = fem.Constant(msh, 1.0e2)   # fibre stiffness

# matrix strain energy (compressible neo-Hookean model)
psi_m = (mu / 2.0) * (I1 - 3.0) - mu * ufl.ln(J) + (lmbda / 2.0) * (ufl.ln(J))**2

# fibre strain energy
psi_f = c4f * (L4f - 1.0)**2

# fibre-matrix interaction strain energy
psi_mf = b4Ff * (J4f**2 - 1.0)

# 1st Piola Kirchhoff-like stress tensors
Pm = ufl.diff(psi_m, F)
Pmf = ufl.diff(psi_mf, F)

# reference 1st Piola-Kirchhoff stress
P0mf = ufl.replace(ufl.algorithms.apply_derivatives.apply_derivatives(ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(Pmf)), {af:Af, F:ufl.Identity(d)})

# ensure Pmf = 0 in reference configuration
Pmf = Pmf - P0mf

# director stress vectors
df = ufl.diff(psi_f, af)
dmf = ufl.diff(psi_mf, af)

# reference director stress
d0mf = ufl.replace(ufl.algorithms.apply_derivatives.apply_derivatives(ufl.algorithms.apply_algebra_lowering.apply_algebra_lowering(dmf)), {af:Af})

# ensure dmf = 0 in reference configuration
dmf = dmf - d0mf

########################################################################################################################
# Variational principle and solver initiatisation

# Variational formulation
metadata = {"quadrature_degree": 4}
dx = ufl.Measure("dx", domain=msh, metadata=metadata)

# Define form Pi (we want to find u such that Pi(u) = 0)
Pi = ufl.inner(Pm, ufl.grad(v)) * dx + ufl.inner(Pmf, ufl.grad(v)) * dx + ufl.inner(df, yf) * dx + ufl.inner(dmf, yf) * dx

########################################################################################################################
# Initialize nonlinear problem and solver
problem = fem.petsc.NonlinearProblem(Pi, uwf, bcs)

# Initialise Newton Rhapson solver
solver = nls.petsc.NewtonSolver(msh.comm, problem)

# solved primary variables
u_h = uwf.sub(0)
wf_h = uwf.sub(1)


########################################################################################################################
# Incremental application of the traction loading
for n in range(1, 10):
    incremented_x_displacement_expr.t = n
    num_its, converged = solver.solve(uwf)
1 Like