Hello!

As the title says, I get some trouble with the splitting function lhs() and rhs() when I try to split the followed variational formulation :

```
f_u = inner( sigma(u,g , G, K0,ndim) , epsilon(v) ) * dx + dot( tn, v ) * ds
a_u = lhs(f_u)
L_u = rhs(f_u)
```

with function define here :

```
# degradation function
def gk( alpha,k = 1e-6):
return (1 - alpha )**2 + k
#shape function
def w(alpha):
return alpha**2
def epsilon(u):
return sym(nabla_grad(u))
def ElasticEnergyDensityPlus(u, K0, G):
return 0.5* K0*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2 + G *inner(dev(epsilon(u)),dev(epsilon(u)))
def sigma(u ,g , G, K0,ndim):
"""
Ko = initial stiffness
G = shear modulus
g = degradation function = gk()
u = displacement
"""
Kd = conditional(gt(tr(epsilon(u)),0),g*K0, K0)
return (2. * Kd * tr(epsilon(u))*Identity(ndim) + G *g * dev(epsilon(u)))
```

Through the stress expression, I try to define the compression/traction distinction of the material behavior. In function of the sign of tr(epsilon), the stiffness coefficient Kd is weigted by a degradation function g

When I run, I get this error : Found Argument in , this is an invalid expression.

Traceback (most recent call last):

File â€śâ€¦/Main_SinglNotc_E2.pyâ€ť, line 216, in

a_u = lhs(f_u)

This is the beginning of my code :

```
from fenics import *
import numpy as np
from Function import *
import matplotlib.pyplot as plt
# ------------------
# Parameters
# ------------------
ERROR = 40
set_log_level(ERROR) # log level
parameters.parse() # read paramaters from command line
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs" # quadrature is deprecated
"""
=========================================================
Time Discretization
=========================================================
"""
T = 1.0 # final time
num_steps = 1000 # number of time steps
dt = T / num_steps # time step size
time_step = np.linspace(0, 1.0, 300)
"""
=========================================================
Create mesh
=========================================================
"""
nx = 50
ny = 50
mesh = UnitSquareMesh(nx, ny,'crossed')
"""
=========================================================
Space discretization
=========================================================
"""
degreeU = 1
degreeAlpha = 1
degreeH = 0
V = VectorFunctionSpace(mesh,'Lagrange', degreeU)
V_alpha = FunctionSpace(mesh,'Lagrange', degreeAlpha)
V_H = FunctionSpace(mesh,'DG', degreeH)
Tens = TensorFunctionSpace(mesh, 'Lagrange', 1)
"""
=========================================================
Define Trial and Test Function
Define functions for solutions at current times steps
=========================================================
"""
u,v,u_ = TrialFunction(V), TestFunction(V), Function(V)
alpha,beta,alpha_ = TrialFunction(V_alpha), TestFunction(V_alpha), Function(V_alpha)
"""
===================================================================
Physical parameters
Parameters
___________________________________________________________________
E= Yound modulus nu : poisson coef
mu=shear modulus K0 : stiffness
mu, lmbda = lamĂ© coeff Gc : critical energy released rate J/mmÂ˛
lw= regularization parameter (internal length) mm
cw=normalized parameter (depends on the shape function w)
for w = phi1Â˛ cw = 2
w = phi1 cw = 3/8
===================================================================
"""
E, nu = Constant(1000), Constant(0.3) #[MPa]
G = Constant(E / (2.0*(1.0 + nu))) # [MPa]
lmbda = Constant (E * nu / ((1.0 + nu)* (1.0 - 2.0*nu)))
K0, Gc = Constant(lmbda+2./3.*G), Constant(0.3)
cw, lw = Constant(2), Constant(0.1) # [mm]
Y_seuil = 0#Constant(Gc/(2*lw))
"""
=========================================================
Boundary condition
=========================================================
"""
#loading : linear displacement (ux,uy) = (0,t/4)
Disp=Expression( ('0.0', 't/4.'),t = 0.0, degree = 1 )
# boolean boundary function
boundary_bottom = 'near(x[1], 0)'
boundary_up = 'near(x[1], 1)'
#Crach_boundary = 'near(x[1], 0.5)' && 'between(x[0], (0., 0.3)) '
class Precrack(SubDomain):
def inside(self, x, on_boundary):
return (between(x[0], (0., 0.3)) and near(x[1], 0.5) ) # and : &&
# Initialize mesh function for exterior domains / corresponding tags to the boundaries
#Up, bottom and crack domain
crack = Precrack()
bottom = CompiledSubDomain(boundary_bottom)
up = CompiledSubDomain(boundary_up)
boundaries = MeshFunction("size_t", mesh,1)
boundaries.set_all(0)
bottom.mark(boundaries, 1) # mark left as 1
up.mark(boundaries, 2) # mark right as 2
crack.mark(boundaries, 3) # mark right as 2
ds = Measure("ds",subdomain_data=boundaries) # left: ds(1), right: ds(2), crack : ds(3), ds(0) all the edges
bcu_bottom = DirichletBC( V , Constant((0.,0.)), boundary_bottom)
bcu_up = DirichletBC( V , Disp, boundary_up )
bcu = [bcu_bottom, bcu_up]
bcp_bottom = DirichletBC(V_alpha, Constant(0.), boundary_bottom)
bcp_up = DirichletBC(V_alpha, Constant(0.), boundary_up)
bcp_crack = DirichletBC(V_alpha, Constant(1.), boundaries,3)
bcp = [bcp_bottom, bcp_up, bcp_crack]
"""
=========================================================
Define Variational problem
=========================================================
"""
# INITIAL CONDITION
u_n = interpolate(Constant((0.,0.)),V)
alpha_n = interpolate(Constant(0.), V_alpha)
H_n = interpolate(Expression("0.0",degree=0), V_H)
H_n_prev = interpolate(Expression("0.0",degree=0), V_H)
#---------------- Mechanical equation
g = gk(alpha_n,1.0e-6) # type : ufl.algebra.Sum
tn = Constant((0., 0.)) # external stress
def sigma(u ,g , G, K0,ndim):
Kd = conditional(gt(tr(epsilon(u)),0),g*K0, K0)
return (2. * Kd * tr(epsilon(u))*Identity(ndim) + G *g * dev(epsilon(u)))
#--------------- Phase Field equation
f_u = inner( sigma(u,g , G, K0,ndim) , epsilon(v) ) * dx + dot( tn, v ) * ds
a_u = lhs(f_u)
L_u = rhs(f_u)
```

I know that the error comes from the conditional expression, but I donâ€™t know what can I do to replace it.

Thanks for your help !