Regarding the problem in User Expression

General
1

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,

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
Screenshot from 2022-11-08 08-39-091464×188 47.6 KB

Thank you
Manish Kumar

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

3

@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:

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
Screenshot from 2022-11-09 18-34-19734×170 33.7 KB

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

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
Screenshot from 2022-11-09 18-40-031445×171 39 KB

Please help me to fix the issue.

Thank you,
Manish

4

Try perhaps

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

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

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

6

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

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
Screenshot from 2022-11-10 08-21-481187×524 72.2 KB

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

Thank you,
Manish

7

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.

8

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

9

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

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.

10

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