# Regarding the problem in User Expression

**URL:** <https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703>\
**Category:** General\
**Created:** [November 8, 2022, 7:41am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703 "2022-11-08T07:41:25Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![manish\_07](https://avatars.discourse-cdn.com/v4/letter/m/e56c9b/32.png) [@manish\_07](https://fenicsproject.discourse.group/u/manish_07)\
**Post date:** [November 8, 2022, 7:41am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/1 "2022-11-08T07:41:25Z")

</div>

Hello everyone,

I am trying to implement the variable stiffness i.e. change the stiffness after some load step in FEniCS.To implement this, I am using User Expression, it was giving me the error. The same expression without User Expression operation works fine. Please help me to fix the issue, I am new to the FEniCS.  
The minimal code is as follow,

```auto
from fenics import *    
from mshr import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np
from ufl import replace
import sys, os, shutil, math
from numpy.linalg import eig
from dolfin import *
from numpy import array
import time
from ufl import eq
from ufl import le
from ufl import ge

# Define Material Properties
E_1 = 22700. # Young's Modulus
nu_1 = 0.2 # Poisson ratio
mu = E_1/(2.0*(1.0 + nu_1)) #lame's constant 
lmbda = (E_1 * nu_1)/((1.0 - 2.0*nu_1)*(1.0 + nu_1)) #lame's constant

# Mesh Data
mesh = Mesh('mesh_Ln5.xml') # Import mesh
ndim = mesh.topology().dim() # Dimension of a geometry

thickness = 50.

# Create function space for 2D elasticity + Damage
V_u = VectorFunctionSpace(mesh, 'CG', 1) # Vector function space for displacement field

# Define the function, test and trial fields (alpha is phase field)
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)

u.rename('displacement','displacement')

tol = 1E-5 # Tolerance for boundary identification
tol_v = 1e-20 # Tolerance for square-root calculation
def boundary_D_l(x, on_boundary): # Find bottom edge of specimen
    return x[0] >= 54.0 - tol and x[0] <= 56.0 - tol and x[1] <= 0.0 + tol

def boundary_D_r(x, on_boundary): # left lower corner to fix
    return x[0] >= 654.0 - tol and x[0] <= 656.0 - tol and x[1] <= 0.0 + tol

def boundary_D_t(x, on_boundary): # Find location of traction/displacement
    return x[1] >= 150.0-tol and x[0] >= 354.98-tol and x[0] <= 355.02+tol 
        
bcl = DirichletBC(V_u, Constant((0.0 ,0.0)), boundary_D_l, method='pointwise') # Apply boundary condition

bcr = DirichletBC(V_u.sub(1), Constant(0.0), boundary_D_r, method='pointwise') # Apply boundary condition

  
# Calculation of positive part of strain
def eps_positive(u):    
    A = sym(grad(u))
    a = A[0,0]
    b = A[0,1]
    c = A[1,0]
    d = A[1,1]
    eig_1 = ((tr(A) + sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    eig_2 = ((tr(A) - sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    phi_1 = (eig_1 - b - d)/(a + c - eig_1)
    phi_2 = (eig_2 - b - d)/(a + c - eig_2)

    eig_v_1 = [phi_1/sqrt(phi_1**2 + 1), 1/sqrt(phi_1**2 + 1)] # Evaluation of Eigen vector
    eig_v_2 = [phi_2/sqrt(phi_2**2 + 1), 1/sqrt(phi_2**2 + 1)] # Evaluation of Eigen vector 

    # Positive Strain Tensor
    sn_P = 0.5*(eig_1 + abs(eig_1))*np.outer(eig_v_1,eig_v_1) + 0.5*(eig_2 + abs(eig_2))*np.outer(eig_v_2,eig_v_2)
    sn_1 = as_matrix(sn_P.tolist())
    return sn_1

# Calculation of negative part of strain
def eps_negative(u):    
    A = sym(grad(u))
    a = A[0,0]
    b = A[0,1]
    c = A[1,0]
    d = A[1,1]
    eig_1 = ((tr(A) + sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    eig_2 = ((tr(A) - sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    phi_1 = (eig_1 - b - d)/(a + c - eig_1)
    phi_2 = (eig_2 - b - d)/(a + c - eig_2)

    eig_v_1 = [phi_1/sqrt(phi_1**2 + 1), 1/sqrt(phi_1**2 + 1)] # Evaluation of Eigen vector
    eig_v_2 = [phi_2/sqrt(phi_2**2 + 1), 1/sqrt(phi_2**2 + 1)] # Evaluation of Eigen vector

         
    # Negative Strain Tensor
    sn_N = 0.5*(eig_1 - abs(eig_1))*np.outer(eig_v_1,eig_v_1) + 0.5*(eig_2 - abs(eig_2))*np.outer(eig_v_2,eig_v_2)
    sn_1 = as_matrix(sn_N.tolist())
    return sn_1

def eps(u):
    """Strain tensor as a function of the displacement"""
    return sym(grad(u))

# expressions_with variable_stiffness
class el_en_1(UserExpression):
    def __init__ (self, u_sol, **kwargs):
        super(). __init__ (**kwargs)
        self.u_sol = u_sol
        
    def eval(self, value, x):
        value = 0.5*lmbda*(0.5*(tr(eps(u_sol)) + abs(tr(eps(u_sol)))))**2 + mu*tr(eps_positive(u_sol)*eps_positive(u_sol))
    def value_shape(self):
        return (2,)

class el_en_2(UserExpression):
    def __init__ (self, u_sol, **kwargs):
        super(). __init__ (**kwargs)
        self.u_sol = u_sol
        
    def eval(self, value, x):
        value = 0.5*lmbda*(0.5*(tr(eps(u_sol)) - abs(tr(eps(u_sol)))))**2 + mu*tr(eps_negative(u_sol)*eps_negative(u_sol))
    def value_shape(self):
        return (2,)
            
ele_ene_1 = el_en_1(u_sol=u,element = V_u.ufl_element())
ele_ene_2 = el_en_2(u_sol=u,element = V_u.ufl_element())
#print(ele_ene_1.shape)
elastic_energy_1 = ele_ene_1#0.5*lmbda*(0.5*(tr(eps(u)) + abs(tr(eps(u)))))**2 + mu*tr(eps_positive(u)*eps_positive(u)) # Positive part of strain energy density
elastic_energy_2 = ele_ene_2#0.5*lmbda*(0.5*(tr(eps(u)) - abs(tr(eps(u)))))**2 + mu*tr(eps_negative(u)*eps_negative(u)) # negative part of strain energy density
elastic_energy = (elastic_energy_1 + elastic_energy_2)*thickness*dx # total strain energy density (bulk energy density)
total_energy = elastic_energy

# First and second directional derivative wrt u
E_u = derivative(total_energy,u,v) # first directional derivatives of total energy with respect of displacement test function v
Jd = derivative(E_u, u, du)

# define loading steps
num_steps = 200 # total number of load steps
total_disp = -0.2 # total applied displacement
steps_loading = 200
disp_steps = total_disp/(steps_loading)

u_R = Expression(('disp_app'),disp_app = disp_steps, n=0., degree=0) # Define loading as expression so that it can be updated for next step
bct = DirichletBC(V_u.sub(1), u_R, boundary_D_t, method='pointwise') # define boundary condition
bc_disp = [bcl, bcr, bct] # Apply boundary conditions

# Solver for displacement field
problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
prm = solver_u.parameters

prm["newton_solver"]["relative_tolerance"] = 5E-1
prm["newton_solver"]["absolute_tolerance"] = 5E-3
prm["newton_solver"]["convergence_criterion"] = "residual"
prm["newton_solver"]["error_on_nonconvergence"] = True
prm["newton_solver"]["linear_solver"] = 'mumps'
prm["newton_solver"]["lu_solver"]["symmetric"] = False 
prm["newton_solver"]["maximum_iterations"] = 2000
prm["newton_solver"]["relaxation_parameter"] = 1.0

prm["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True 

for n in range(num_steps):
    u_R.disp_app = disp_steps*(n+1)
    solver_u.solve()
    

```

The error is as follow,

 ![Screenshot from 2022-11-08 08-39-09](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/3/3ad839bce9ddeda072705f70ac6f3a036c036fe8.png)

Thank you  
Manish Kumar

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 8, 2022, 7:48am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/2 "2022-11-08T07:48:32Z")

</div>

> [@manish\_07](#):
>
> ```auto
> class el_en_1(UserExpression):
> def __init__ (self, u_sol, **kwargs):
> super(). __init__ (**kwargs)
> self.u_sol = u_sol
>         
> def eval(self, value, x):
> value = 0.5*lmbda*(0.5*(tr(eps(u_sol)) + abs(tr(eps(u_sol)))))**2 + mu*tr(eps_positive(u_sol)*eps_positive(u_sol))
> def value_shape(self):
> return (2,)
> 
> class el_en_2(UserExpression):
> def __init__ (self, u_sol, **kwargs):
> super(). __init__ (**kwargs)
> self.u_sol = u_sol
>         
> def eval(self, value, x):
> value = 0.5*lmbda*(0.5*(tr(eps(u_sol)) - abs(tr(eps(u_sol)))))**2 + mu*tr(eps_negative(u_sol)*eps_negative(u_sol))
> def value_shape(self):
> return (2,)
> 
> ```

You cannot use ufl operators such as tr, eps (sym(grad())) etc inside a user expression. You would have to project these expressions into a suitable function space

---

<div class="post-metadata">

**Author:** ![manish\_07](https://avatars.discourse-cdn.com/v4/letter/m/e56c9b/32.png) [@manish\_07](https://fenicsproject.discourse.group/u/manish_07)\
**Post date:** [November 9, 2022, 5:42pm UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/3 "2022-11-09T17:42:20Z")

</div>

@dokken Thank you for your quick response. I tried as you suggested by projecting the ufl operator and then take those as input arguments for the user expression. But it is throwing the error:

_This integral is missing an integration domain._

I think it is due to some problem with functionspace but I am not able to find the problem. Can you please have a look on the code and help me to fix the issue.

Minimal code is as follows:

```auto
from fenics import *    
from mshr import *
import numpy as np
from ufl import replace
from dolfin import *
import time

# Define Material Properties
E_1 = 22700. # Young's Modulus
nu_1 = 0.2 # Poisson ratio
mu = E_1/(2.0*(1.0 + nu_1)) #lame's constant 
lmbda = (E_1 * nu_1)/((1.0 - 2.0*nu_1)*(1.0 + nu_1)) #lame's constant

# Mesh Data
mesh = Mesh('mesh_Ln5.xml') # Import mesh
ndim = mesh.topology().dim() # Dimension of a geometry

thickness = 50.

# Create function space for 2D elasticity + Damage
V_u = VectorFunctionSpace(mesh, 'CG', 1) # Vector function space for displacement field
V_T = TensorFunctionSpace(mesh, 'DG', 0)
V_S = FunctionSpace(mesh, 'DG', 0)

# Define the function, test and trial fields (alpha is phase field)
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)

u.rename('displacement','displacement')

tol = 1E-5 # Tolerance for boundary identification
tol_v = 1e-20 # Tolerance for square-root calculation
def boundary_D_l(x, on_boundary): # Find bottom edge of specimen
    return x[0] >= 54.0 - tol and x[0] <= 56.0 - tol and x[1] <= 0.0 + tol

def boundary_D_r(x, on_boundary): # left lower corner to fix
    return x[0] >= 654.0 - tol and x[0] <= 656.0 - tol and x[1] <= 0.0 + tol

def boundary_D_t(x, on_boundary): # Find location of traction/displacement
    return x[1] >= 150.0-tol and x[0] >= 354.98-tol and x[0] <= 355.02+tol 
        
bcl = DirichletBC(V_u, Constant((0.0 ,0.0)), boundary_D_l, method='pointwise') # Apply boundary condition

bcr = DirichletBC(V_u.sub(1), Constant(0.0), boundary_D_r, method='pointwise') # Apply boundary condition

  
# Calculation of positive part of strain
def eps_positive(u):    
    A = sym(grad(u))
    a = A[0,0]
    b = A[0,1]
    c = A[1,0]
    d = A[1,1]
    eig_1 = ((tr(A) + sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    eig_2 = ((tr(A) - sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    phi_1 = (eig_1 - b - d)/(a + c - eig_1)
    phi_2 = (eig_2 - b - d)/(a + c - eig_2)

    eig_v_1 = [phi_1/sqrt(phi_1**2 + 1), 1/sqrt(phi_1**2 + 1)] # Evaluation of Eigen vector
    eig_v_2 = [phi_2/sqrt(phi_2**2 + 1), 1/sqrt(phi_2**2 + 1)] # Evaluation of Eigen vector 

    # Positive Strain Tensor
    sn_P = 0.5*(eig_1 + abs(eig_1))*np.outer(eig_v_1,eig_v_1) + 0.5*(eig_2 + abs(eig_2))*np.outer(eig_v_2,eig_v_2)
    sn_1 = as_matrix(sn_P.tolist())
    return sn_1

def eps(u):
    """Strain tensor as a function of the displacement"""
    return sym(grad(u))

# expressions_with variable_stiffness
eps_ = project(eps(u),V_T)
eps_p = project(eps_positive(u)*eps_positive(u),V_T)

tr_eps_ = project((tr(eps_) + abs(tr(eps_))),V_S)      
tr_eps_p = project(tr(eps_p),V_S)

class el_en_1(UserExpression):
    def __init__ (self, tr_eps_, tr_eps_p, lmbda, mu, **kwargs):
        super(). __init__ (**kwargs)
        self.tr_eps_sol = tr_eps_
        self.tr_eps_p_sol = tr_eps_p
        self.lmbda_sol = lmbda
        self.mu_sol = mu
        
    def eval(self, value, x):
        value = 0.5*self.lmbda_sol*(0.5*(self.tr_eps_sol))**2 + self.mu_sol*self.tr_eps_p_sol
    def value_shape(self):
        return ()

ele_ene_1 = el_en_1(tr_eps_=tr_eps_, tr_eps_p=tr_eps_p, lmbda = lmbda, mu = mu)

elastic_energy_1 = ele_ene_1
elastic_energy = (elastic_energy_1)*thickness*dx  
total_energy = elastic_energy

# First and second directional derivative wrt u
E_u = derivative(total_energy,u,v) # first directional derivatives of total energy with respect of displacement test function v
Jd = derivative(E_u, u, du)

# define loading steps
num_steps = 200 # total number of load steps
total_disp = -0.2 # total applied displacement
steps_loading = 200
disp_steps = total_disp/(steps_loading)

u_R = Expression(('disp_app'),disp_app = disp_steps, n=0., degree=0) # Define loading as expression so that it can be updated for next step
bct = DirichletBC(V_u.sub(1), u_R, boundary_D_t, method='pointwise') # define boundary condition
bc_disp = [bcl, bcr, bct] # Apply boundary conditions

# Solver for displacement field
problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
prm = solver_u.parameters

prm["newton_solver"]["relative_tolerance"] = 5E-1
prm["newton_solver"]["absolute_tolerance"] = 5E-3
prm["newton_solver"]["convergence_criterion"] = "residual"
prm["newton_solver"]["error_on_nonconvergence"] = True
prm["newton_solver"]["linear_solver"] = 'mumps'
prm["newton_solver"]["lu_solver"]["symmetric"] = False 
prm["newton_solver"]["maximum_iterations"] = 2000
prm["newton_solver"]["relaxation_parameter"] = 1.0

prm["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True 

for n in range(num_steps):
    u_R.disp_app = disp_steps*(n+1)
    solver_u.solve()
    

```

The error is as follows:

 ![Screenshot from 2022-11-09 18-34-19](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/a70989b37a68f3a8a02002c21fc1eae27aff1a05.png)

I tried by providing the element argument also at the time of calling user expression as follows:

```auto
ele_ene_1 = el_en_1(tr_eps_=tr_eps_, tr_eps_p=tr_eps_p, lmbda = lmbda, mu = mu,element = V_u.ufl_element())

```

It throws the error that integrand is tensor as follows:

 ![Screenshot from 2022-11-09 18-40-03](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/2/22e5af55d55c6a03118cebbcf99423a3fd58dc6e.png)

Please help me to fix the issue.

Thank you,  
Manish

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [November 9, 2022, 6:00pm UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/4 "2022-11-09T18:00:18Z")

</div>

Try perhaps

```auto
ele_ene_1 = el_en_1(tr_eps_=tr_eps_, tr_eps_p=tr_eps_p, lmbda = lmbda, mu = mu, domain=mesh)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 9, 2022, 7:05pm UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/5 "2022-11-09T19:05:12Z")

</div>

> [@manish\_07](#):
>
> @dokken Thank you for your quick response. I tried as you suggested by projecting the ufl operator and then take those as input arguments for the user expression. But it is throwing the error:

You need to project the whole expression ,ie. all of  
` value = 0.5*lmbda*(0.5*(tr(eps(u_sol)) - abs(tr(eps(u_sol)))))**2 + mu*tr(eps_negative(u_sol)*eps_negative(u_sol))`  
This means

```python
value = 0.5*lmbda*(0.5*(tr(eps(u_sol)) - abs(tr(eps(u_sol)))))**2 + mu*tr(eps_negative(u_sol)*eps_negative(u_sol))
expr_function = project(value, Compatible_function_space)

```

You cannot do linear operations (as in multiply functions from the same function space together and end up in the same function space).

Also, the integral measure `dx` could take `dx(domain=mesh)` as @nate suggested

---

<div class="post-metadata">

**Author:** ![manish\_07](https://avatars.discourse-cdn.com/v4/letter/m/e56c9b/32.png) [@manish\_07](https://fenicsproject.discourse.group/u/manish_07)\
**Post date:** [November 10, 2022, 7:31am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/6 "2022-11-10T07:31:50Z")

</div>

@dokken Thank you for your suggestions. I tried your suggestions and the updated code is as follows

```auto
from fenics import *    
from mshr import *
import numpy as np
from ufl import replace
from dolfin import *
import time

# Define Material Properties
E_1 = 22700. # Young's Modulus
nu_1 = 0.2 # Poisson ratio
mu = E_1/(2.0*(1.0 + nu_1)) #lame's constant 
lmbda = (E_1 * nu_1)/((1.0 - 2.0*nu_1)*(1.0 + nu_1)) #lame's constant

# Mesh Data
mesh = Mesh('mesh_Ln5.xml') # Import mesh
ndim = mesh.topology().dim() # Dimension of a geometry

thickness = 50.

# Create function space for 2D elasticity + Damage
V_u = VectorFunctionSpace(mesh, 'CG', 1) # Vector function space for displacement field
V_T = TensorFunctionSpace(mesh, 'DG', 0)
V_S = FunctionSpace(mesh, 'DG', 0)

# Define the function, test and trial fields (alpha is phase field)
u, du, v = Function(V_u), TrialFunction(V_u), TestFunction(V_u)

u.rename('displacement','displacement')

tol = 1E-5 # Tolerance for boundary identification
tol_v = 1e-20 # Tolerance for square-root calculation
def boundary_D_l(x, on_boundary): # Find bottom edge of specimen
    return x[0] >= 54.0 - tol and x[0] <= 56.0 - tol and x[1] <= 0.0 + tol

def boundary_D_r(x, on_boundary): # left lower corner to fix
    return x[0] >= 654.0 - tol and x[0] <= 656.0 - tol and x[1] <= 0.0 + tol

def boundary_D_t(x, on_boundary): # Find location of traction/displacement
    return x[1] >= 150.0-tol and x[0] >= 354.98-tol and x[0] <= 355.02+tol 
        
bcl = DirichletBC(V_u, Constant((0.0 ,0.0)), boundary_D_l, method='pointwise') # Apply boundary condition

bcr = DirichletBC(V_u.sub(1), Constant(0.0), boundary_D_r, method='pointwise') # Apply boundary condition

  
# Calculation of positive part of strain
def eps_positive(u):    
    A = sym(grad(u))
    a = A[0,0]
    b = A[0,1]
    c = A[1,0]
    d = A[1,1]
    eig_1 = ((tr(A) + sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    eig_2 = ((tr(A) - sqrt(tr(A)**2-4*det(A) + tol_v))/2) # Evaluation of Eigen value 
    phi_1 = (eig_1 - b - d)/(a + c - eig_1)
    phi_2 = (eig_2 - b - d)/(a + c - eig_2)

    eig_v_1 = [phi_1/sqrt(phi_1**2 + 1), 1/sqrt(phi_1**2 + 1)] # Evaluation of Eigen vector
    eig_v_2 = [phi_2/sqrt(phi_2**2 + 1), 1/sqrt(phi_2**2 + 1)] # Evaluation of Eigen vector 

    # Positive Strain Tensor
    sn_P = 0.5*(eig_1 + abs(eig_1))*np.outer(eig_v_1,eig_v_1) + 0.5*(eig_2 + abs(eig_2))*np.outer(eig_v_2,eig_v_2)
    sn_1 = as_matrix(sn_P.tolist())
    return sn_1

def eps(u):
    """Strain tensor as a function of the displacement"""
    return sym(grad(u))

# expressions_with variable_stiffness
eps_ = project(eps(u),V_T)
eps_p = project(eps_positive(u)*eps_positive(u),V_T)

tr_eps_ = project((tr(eps_) + abs(tr(eps_))),V_S)      
tr_eps_p = project(tr(eps_p),V_S)

class el_en_1(UserExpression):
    def __init__ (self, tr_eps_, tr_eps_p, lmbda, mu,domain, **kwargs):
        super(). __init__ (**kwargs)
        self.tr_eps_sol = tr_eps_
        self.tr_eps_p_sol = tr_eps_p
        self.lmbda_sol = lmbda
        self.mu_sol = mu
        self.domain = domain
    def eval(self, value, x):
        value = 0.5*self.lmbda_sol*(0.5*(self.tr_eps_sol))**2 + self.mu_sol*self.tr_eps_p_sol
    def value_shape(self):
        return ()

ele_ene_1 = el_en_1(tr_eps_=tr_eps_, tr_eps_p=tr_eps_p, lmbda = lmbda, mu = mu,domain=mesh)

elastic_energy_1 = ele_ene_1
elastic_energy = (elastic_energy_1)*thickness*dx(domain=mesh)
total_energy = elastic_energy

# First and second directional derivative wrt u
E_u = derivative(total_energy,u,v) # first directional derivatives of total energy with respect of displacement test function v
Jd = derivative(E_u, u, du)

# define loading steps
num_steps = 200 # total number of load steps
total_disp = -0.2 # total applied displacement
steps_loading = 200
disp_steps = total_disp/(steps_loading)

u_R = Expression(('disp_app'),disp_app = disp_steps, n=0., degree=0) # Define loading as expression so that it can be updated for next step
bct = DirichletBC(V_u.sub(1), u_R, boundary_D_t, method='pointwise') # define boundary condition
bc_disp = [bcl, bcr, bct] # Apply boundary conditions

# Solver for displacement field
problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
prm = solver_u.parameters

prm["newton_solver"]["relative_tolerance"] = 5E-1
prm["newton_solver"]["absolute_tolerance"] = 5E-3
prm["newton_solver"]["convergence_criterion"] = "residual"
prm["newton_solver"]["error_on_nonconvergence"] = True
prm["newton_solver"]["linear_solver"] = 'mumps'
prm["newton_solver"]["lu_solver"]["symmetric"] = False 
prm["newton_solver"]["maximum_iterations"] = 2000
prm["newton_solver"]["relaxation_parameter"] = 1.0

prm["newton_solver"]["krylov_solver"]["nonzero_initial_guess"] = True 

for n in range(num_steps):
    u_R.disp_app = disp_steps*(n+1)
    solver_u.solve()
    

```

but it is throwing a new error as follows:

 ![Screenshot from 2022-11-10 08-21-48](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/1/1112b387fb7db97a8947876de3061422a3b2c51d.png)

The residual error after first iteration is different than the approach without user expression.

Thank you,  
Manish

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 10, 2022, 7:46am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/7 "2022-11-10T07:46:21Z")

</div>

Wait a second. You are using this as input to the residual F. That cannot work. Why do you want this to be a UserExpression at all?

I dont see why this cannot be implemented directly with ufl expressions (no user expression).

There is a lot going on in your code, and I dont have time to parse it all. Please try to make a minimal example of what you want to do. It does not need to involve a solve.

---

<div class="post-metadata">

**Author:** ![manish\_07](https://avatars.discourse-cdn.com/v4/letter/m/e56c9b/32.png) [@manish\_07](https://fenicsproject.discourse.group/u/manish_07)\
**Post date:** [November 10, 2022, 8:20am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/8 "2022-11-10T08:20:59Z")

</div>

I want to model loading and unloading. During unloading, the stiffness is a function of a solution parameter. So, my idea was to use the User Expression and change the stiffness during unloading.

Thank you,  
Manish

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 10, 2022, 10:09am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/9 "2022-11-10T10:09:50Z")

</div>

But why cant you write that directly in your variational form without an expression, i.e.

```python
elastic_energy_1 = 0.5*lmbda*(0.5*(tr(eps(u)) + abs(tr(eps(u)))))**2 + mu*tr(eps_positive(u)*eps_positive(u))
elastic_energy_2 = 0.5*lmbda*(0.5*(tr(eps(u)) - abs(tr(eps(u)))))**2 + mu*tr(eps_negative(u)*eps_negative(u)) # negative part of strain energy density
elastic_energy = (elastic_energy_1 + elastic_energy_2)*thickness*dx # total strain energy density (bulk energy density)
total_energy = elastic_energy

```

etc.  
As `u` is the unknown you are solving for, this should incorporate it in the varitational formulation.

---

<div class="post-metadata">

**Author:** ![manish\_07](https://avatars.discourse-cdn.com/v4/letter/m/e56c9b/32.png) [@manish\_07](https://fenicsproject.discourse.group/u/manish_07)\
**Post date:** [November 10, 2022, 10:32am UTC](https://fenicsproject.discourse.group/t/regarding-the-problem-in-user-expression/9703/10 "2022-11-10T10:32:25Z")

</div>

Thank you for your suggestions. I will work in this way.