Mixed element method in Fenics for nonlinear fracture problem

General
1

Greetings, I am trying to run a nonlinear fracture simulation using the mixed element method in fenics and no matter how many corrects I implemented using the advice from this forum, I keep getting the following error:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range).
*** Where:   This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1617882212586/work/dolfin/dolfin/la/PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  7f46aeb0b296da5bbb1fb0845822a72ab9b09c55
*** -------------------------------------------------------------------------

Would really appreciate it if someone knew what the problem is. I attached the code below thank you!

# Janel Chua
# Able to handle supersonic (with respect to pressure wave speed)
# in bridges: cd /ocean/projects/dmr120014p/songlinc/
#########################################################################################
# Preliminaries and mesh
from dolfin import *
import numpy as np
# mesh_half = Mesh('Trial14_coarseEdges8.xml')
# mesh_half = Mesh('Trial31.xml')
mesh_half = Mesh('Trial17.xml')

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Manually introduce the material parameters
class GC(UserExpression):
    def set_Gc_values(self, Gc_0, Gc_1):
        self.Gc_0, self.Gc_1 = Gc_0, Gc_1
    def eval(self, value, x):
        "Set value[0] to value at point x"
        tol = 1E-14
        if x[1] >= 0.015 + tol:
            value[0] = self.Gc_0
        elif x[1] <= -0.015 + tol:
            value[0] = self.Gc_0
        else:
            value[0] = self.Gc_1 #middle layer
# Initialize Gc
Gc = GC()
Gc.set_Gc_values(32, 0.1) #N/mm 

# STEEL
l     = 0.015#*10**(-3) #mm
E     = 206*1000 #MPa = MKg.m^-1.s^-2 #we are using N/mm^2
nu    = 0.3 #Poisson's ratio
mu    = Constant(E / (2.0*(1.0 + nu))) #80.7692e3# MPa #we are using N/mm^2
lmbda = Constant(E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))) #121.1538e3# MPa #we are using N/mm^2
mu2   = 1e-10 #this mu is used for the evolution equation
# Mass density
rho   = Constant(7.8*10**(-9)) # Mkg.m^-3  #we are using (N.s^2)/mm^4
#----------------------------------------------------------------------#

# Generalized-alpha method parameters
alpha_m = Constant(0.2)
alpha_f = Constant(0.4)
gamma   = Constant(0.5+alpha_f-alpha_m)
beta    = Constant((gamma+0.5)**2/4.)

# Time-stepping parameters
T       = 20 
Nsteps  = 100 
dt = Constant(T/Nsteps)

tract = Constant((0, 5)) # Vertical traction

#########################################################################################
# Define Function Spaces and mixed elements
FE_A  = VectorElement("CG", mesh_half.ufl_cell(), 1)
FE_U  = VectorElement("CG", mesh_half.ufl_cell(), 1) # 1 originally # displacement
FE_Phi  = FiniteElement("CG", mesh_half.ufl_cell(), 1) # phi

Z = FunctionSpace(mesh_half, MixedElement([FE_A, FE_U, FE_Phi]))
dU = TrialFunction(Z)
(A, U, Phi) = split(dU)
U_ = TestFunction(Z)
(V, W, Q) = split(U_)

old = Function(Z)
(Aold, Uold, Phiold) = split(old)
new = Function(Z)

strain = TensorFunctionSpace(mesh_half, "CG", 1)
strain_total = Function(strain, name='Strain')
sigma_fs = TensorFunctionSpace(mesh_half, "CG", 1)
stress_elas = Function(sigma_fs, name='Stress')
VV = FunctionSpace(mesh_half, 'CG', 1)
normSquared = Function(VV)

# Defining velocity of moving defect
vel = as_vector([3000000,0]) # mm/s
velx = vel[0] 

# Rayleigh damping coefficients
eta_k = Function(VV)
class etakValues(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        etakValue = Constant(1*10**(-9)) # etak's maximum value 
        etakY = 0.2 # In this case, y position where you would like the value to be turned on. Applies to Pos and Neg Y
        etakPosX = 3
        etakNegX = 0 #-2.5 # if you want it to be in the positive region you have to add a -
        l2 = 0.1
        baseline = Constant(0) # Constant(1*10**(-9)) # this is the value in the middle section

        def H(u):
            return 0.5*( 1 + ((np.exp(2*((u)/l2)) - 1)/(np.exp(2*((u)/l2))+1)) )
        # value[0] = (-etakValue)*H(x[0] + etakNegX)*H(x[1] + etakY)*(1-H(x[0] - etakPosX))*(1-H(x[1] - etakY))\
        # + etakValue + baseline # Good
        # value[0] = baseline 
        value[0] = etakValue
eta_k.interpolate(etakValues())

# Phi initialisation
class PhiInitialCondition(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        if abs(x[1]) < 1e-02 and x[0] <= 0.1:
            value[0] = 1.0
        else:
            value[0] = 0.0
Phiold = interpolate(PhiInitialCondition(), Z.sub(2).collapse())

#########################################################################################
# Boundary conditions
def top(x,on_boundary):
    return near(x[1],1) and on_boundary #(x[0], 1)
def bot(x,on_boundary):
    return near(x[1],-1) and on_boundary #(x[0], -1)
def left(x,on_boundary):
    return near(x[0],-0.5) and on_boundary #(-0.5, x[1])
def right(x,on_boundary):
    return near(x[0],5) and on_boundary #(5, x[1])

bcright = DirichletBC(Z.sub(1).sub(0), Constant(0), right) #Right displacement loaded: u_1
bcright2 = DirichletBC(Z.sub(1).sub(1), Constant(0),  right) #Right displacement loaded: u_2
phiright = DirichletBC(Z.sub(2), Constant(0), right)

# Dirichlet Boundary Condition used:
bc_U = [bcright, bcright2, phiright]  

n = FacetNormal(mesh_half)
# Create mesh function over the cell facets
boundary_subdomains = MeshFunction("size_t", mesh_half, mesh_half.topology().dim() - 1)
boundary_subdomains.set_all(0)

AutoSubDomain(top).mark(boundary_subdomains, 1) # top boundary
AutoSubDomain(bot).mark(boundary_subdomains, 2) # bottom boundary
AutoSubDomain(left).mark(boundary_subdomains, 3) # left boundary
AutoSubDomain(leftTopHalf).mark(boundary_subdomains, 31) # leftTop boundary
AutoSubDomain(leftBotHalf).mark(boundary_subdomains, 32) # leftBot boundary
AutoSubDomain(right).mark(boundary_subdomains, 4) # right boundary

# Define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)

#########################################################################################
# Constituive functions
def epsilon(U):
    return sym(grad(U))
def epsilonDot(a):
    return -sym(grad(a))
def sigma(U, a, Phi):
    eta_e = 1e-10
    elastic = ((1.0 - Phi)**2 + eta_e)*(2.0*mu*epsilon(U)+lmbda*tr(epsilon(U))*Identity(len(U)))
    dissipative = ((1.0 - Phi)**2 + eta_e)*(eta_k*( 2.0*mu*epsilonDot(a)+lmbda*tr(epsilonDot(a))*Identity(len(U)) ))
    return (elastic + dissipative)
def psi_elastic(U):
    return 0.5*lmbda*(tr(epsilon(U)))**2 + mu*inner(epsilon(U),epsilon(U)) # isotropic linear elastic


# Mass term weak form
def m(U, a, W, vel):
    return rho*inner(a, div(outer(W, vel)) )*dx
# Elastic stiffness weak form
def k(U, a, Phi, W):
    return inner(sigma(U, a, Phi), grad(W))*dx 
# Boundary Terms
def boundaryLeft(U, W, velx, a, Phi):
    stress = sigma(U, a, Phi)
    e1 = as_vector([1,0])
    e2 = as_vector([0,1])
    n = as_vector([-1,0])
    return -rho*velx*a[0]*dot(W,e1)*n[0] - rho*velx*a[1]*dot(W,e2)*n[0] 
def boundary22(Phi, W, tract): # For traction boundary
    eta_e = 1e-10
    return ((1.0 - Phi)**2 + eta_e)*dot(W, tract) 

#########################################################################################
# Weak form for a -------------------------------------------------------------#
E_A = ( inner(dot(grad(Uold),vel), V) - inner(Aold, V) )*dx

# Weak form for u -------------------------------------------------------------#
# Nonlinear
E_U = m(Uold, Aold, W, vel) \
      - k(Uold, Aold, Phiold, W ) \
      + ( boundaryLeft(Uold, W, velx, Aold, Phiold ) )*dss(3) \
      + ( boundary22(Phiold, W, tract) )*dss(31) \
      + ( boundary22(Phiold, W, -tract) )*dss(32) 

# Weak form for Phi... For the fracture case ----------------------------------#
E_Phi = (-mu2*inner(grad(Phiold),vel)*Q + Gc*l*inner(grad(Phiold),grad(Q)) + ((Gc/l)+2.0*psi_elastic(Uold))*inner(Phiold,Q) - 2.0*psi_elastic(Uold)*Q)*dx

form = E_A + E_U + E_Phi

# Nonlinear Solver
J = derivative(form, old, dU)
# Solve problem using Newton's method
problem = NonlinearVariationalProblem(form, old, bc_U, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
info(prm, True)
prm["nonlinear_solver"] = "newton"
prm["newton_solver"]["absolute_tolerance"] = 1E-8
prm["newton_solver"]["relative_tolerance"] = 1E-8
prm["newton_solver"]["maximum_iterations"] = 100

#########################################################################################
# Initialization of the iterative procedure and output requests
time = np.linspace(0, T, Nsteps+1)
tol = 1e10 

#Name of file used for storage
store_Phi = File ("mydata/Phi.pvd")
store_u = File ("mydata/u.pvd")
store_a = File ("mydata/a.pvd")
store_stress_elas = File ("mydata/stress_elas.pvd")

#########################################################################################
# Looping through time here.
for (i, dt) in enumerate(np.diff(time)):

    t = time[i+1]
    print("Time: ", t)
    iter = 0
    err = 1e11

    while err > tol:
        iter += 1
        # Solve for new displacement  
        # Nonlinear Solver
        solver.solve()
        # Split the solution and plot results
        (Anew, Unew, Phinew) = old.split()

        # Calculate error
        err_a = errornorm(Anew,Aold,norm_type = 'l2',mesh = None)
        err_u = errornorm(Unew,Uold,norm_type = 'l2',mesh = None)
        err_Phi = errornorm(Phinew,Phiold,norm_type = 'l2',mesh = None)
        err = max(err_a, err_u, err_Phi)
        print(err_a, err_u, err_Phi)

        # Update new fields in same timestep with new calculated quantities
        Aold.assign(Anew)
        Uold.assign(Unew)
        Phiold.assign(Phinew)

        print ('Iterations:', iter, ', Total time', t, ', error', err)
	
        if err < tol:
            # Update old fields from previous timestep with new quantities
            Aold.assign(Anew)
            Uold.assign(Unew)
            Phiold.assign(Phinew)

            print ('err<tol :D','Iterations:', iter, ', Total time', t, ', error', err)

            if round(t*1e1) % 2 == 0: # each saved data point is 2e-9s
                store_a << Aold
                store_u << Uold #mm
                store_Phi << Phiold

                stress_elas.assign(project( sigma(Uold, Aold, Phiold) ,sigma_fs, solver_type="cg", preconditioner_type="amg")) # 1MPa = 1N/mm^2
                store_stress_elas << stress_elas

                File('mydata/saved_mesh.xml') << mesh_half
                File('mydata/saved_a.xml') << Aold
                File('mydata/saved_u.xml') << Uold
                print ('Iterations:', iter, ', Total time', t, 'Saving datapoint')
 	    
print ('Simulation completed') 
#########################################################################################
2

Please make a minimal reproducible example.
There is no way we can reproduce the behavior of your code, as, we do not have the mesh.

One thing that stands out is that you overload Phiold from

which could cause the issue.

Secondly, your problem is far from minimal. Please try to reduce the complexity, as there is alot of code here, and takes a while to digest.

Finally,

Please search the forum for this error message as well,
as you find it in many recent threads, such as:

which points to the issue I’ve highlighted above.

3

Hi Dokken,

Thank you so much for getting back to me!
Based on your comments, I commented out the interpolation of Phiold like you suggested and also changed my nonlinear solver to the following:

# Nonlinear Solver
J = derivative(form, old, dU)
# Solve problem using Newton's method
problem = NonlinearVariationalProblem(form, old, bc_U, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
info(prm, True)
prm["nonlinear_solver"] = "newton"
prm["newton_solver"]["linear_solver"] = "mumps"
prm["newton_solver"]["absolute_tolerance"] = 1E4
prm["newton_solver"]["relative_tolerance"] = 1E4
prm["newton_solver"]["maximum_iterations"] = 100

There are no more runtime errors after those two changes! Thank you so much for your help!

4

Hi dokken,

Sorry to bother you again, in my last reply there are no runtime errors, but i found that the dirichlet boundary conditions are not applied the way i coded it to… For instance, in the code, bcleft and phileft simply do not show up when i check for it in paraview…

I made the code much more minimal this time and included a mesh. Would really appreciate any suggestions on why my dirichlet BC are not applied.

# Janel Chua
#########################################################################################
# Preliminaries and mesh
from dolfin import *
import numpy as np
from math import isclose
# Define the corner points and resolution
x0, y0 = -0.5, -1   # Bottom-left corner (start point)
x1, y1 = 4.5, 1.0   # Top-right corner (end point)
nx, ny = 500, 200     # Number of divisions in the x and y directions

# Create the rectangular mesh
mesh_half = RectangleMesh(Point(x0, y0), Point(x1, y1), nx, ny)

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Manually introduce the material parameters
class GC(UserExpression):
    def set_Gc_values(self, Gc_0, Gc_1):
        self.Gc_0, self.Gc_1 = Gc_0, Gc_1
    def eval(self, value, x):
        "Set value[0] to value at point x"
        tol = 1E-14
        if x[1] >= 0.015 + tol:
            value[0] = self.Gc_0
        elif x[1] <= -0.015 + tol:
            value[0] = self.Gc_0
        else:
            value[0] = self.Gc_1 #middle layer
# Initialize Gc
Gc = GC()
Gc.set_Gc_values(32, 0.1) #N/mm 

# STEEL
l     = 0.015#*10**(-3) #mm
E     = 206*1000 #MPa = MKg.m^-1.s^-2 #we are using N/mm^2
nu    = 0.3 #Poisson's ratio
mu    = Constant(E / (2.0*(1.0 + nu))) #80.7692e3# MPa #we are using N/mm^2
lmbda = Constant(E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))) #121.1538e3# MPa #we are using N/mm^2
mu2   = 1e-10 #this mu is used for the evolution equation
# Mass density
rho   = Constant(7.8*10**(-9)) # Mkg.m^-3  #we are using (N.s^2)/mm^4
#----------------------------------------------------------------------#

# Time-stepping parameters
T       = 20 
Nsteps  = 100 
dt = Constant(T/Nsteps)

tract = Constant((0, 5)) # Vertical traction
DOLFIN_EPS = 1.0e-6
#########################################################################################
# Define Function Spaces and mixed elements
FE_A  = VectorElement("CG", mesh_half.ufl_cell(), 1)
FE_U  = VectorElement("CG", mesh_half.ufl_cell(), 1) # 1 originally # displacement
FE_Phi  = FiniteElement("CG", mesh_half.ufl_cell(), 1) # phi

Z = FunctionSpace(mesh_half, MixedElement([FE_A, FE_U, FE_Phi]))
print(Z.num_sub_spaces()) 
dU = TrialFunction(Z)
(A, U, Phi) = split(dU)
U_ = TestFunction(Z)
(V, W, Q) = split(U_)

old = Function(Z)
(Aold, Uold, Phiold) = old.split(True)

sigma_fs = TensorFunctionSpace(mesh_half, "CG", 1)
stress_elas = Function(sigma_fs, name='Stress')
VV = FunctionSpace(mesh_half, 'CG', 1)

# Defining velocity of moving defect
vel = as_vector([3000000,0]) # mm/s
velx = vel[0] 

# Rayleigh damping coefficients
eta_k = Function(VV)
class etakValues(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        etakValue = Constant(1*10**(-9)) # etak's maximum value 
        value[0] = etakValue
eta_k.interpolate(etakValues())

# Phi initialisation
class PhiInitialCondition(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        if abs(x[1]) < 1e-02 and x[0] <= 0.1:
            value[0] = 1.0
        else:
            value[0] = 0.0
Phiold.interpolate(PhiInitialCondition())

#########################################################################################
# Boundary conditions
def left(x,on_boundary):
    return near(x[0],-0.5,DOLFIN_EPS) and on_boundary #(-0.5, x[1])
def right(x,on_boundary):
    return near(x[0],4.5,DOLFIN_EPS) and on_boundary #(4.5, x[1])
def imposedPhi(x,on_boundary):
    return near(x[0],-0.5,DOLFIN_EPS) and on_boundary

bcright = DirichletBC(Z.sub(1).sub(0), Constant(0), right) #Right displacement loaded: u_1
bcright2 = DirichletBC(Z.sub(1).sub(1), Constant(0),  right) #Right displacement loaded: u_2
stretch1 = Expression(("0.1"),degree=1)
bcleft= DirichletBC(Z.sub(1).sub(0), stretch1, left) #Left displacement loaded: u_1
phiright = DirichletBC(Z.sub(2), Constant(0), right)
phileft = DirichletBC(Z.sub(2), Constant(1.0), imposedPhi)

# Dirichlet Boundary Condition used:
bc_U = [bcright, bcright2, bcleft, phiright, phileft]  

n = FacetNormal(mesh_half)

#########################################################################################
# Constituive functions
def epsilon(U):
    return sym(grad(U))
def epsilonDot(a):
    return -sym(grad(a))
def sigma(U, a, Phi):
    eta_e = 1e-10
    elastic = ((1.0 - Phi)**2 + eta_e)*(2.0*mu*epsilon(U)+lmbda*tr(epsilon(U))*Identity(len(U)))
    dissipative = ((1.0 - Phi)**2 + eta_e)*(eta_k*( 2.0*mu*epsilonDot(a)+lmbda*tr(epsilonDot(a))*Identity(len(U)) ))
    return (elastic + dissipative)
def psi_elastic(U):
    return 0.5*lmbda*(tr(epsilon(U)))**2 + mu*inner(epsilon(U),epsilon(U)) # isotropic linear elastic

# Mass term weak form
def m(U, a, W, vel):
    return rho*inner(a, div(outer(W, vel)) )*dx
# Elastic stiffness weak form
def k(U, a, Phi, W):
    return inner(sigma(U, a, Phi), grad(W))*dx 

#########################################################################################
# Weak form for a -------------------------------------------------------------#
E_A = ( inner(dot(grad(Uold),vel), V) - inner(Aold, V) )*dx

# Weak form for u -------------------------------------------------------------#
# Nonlinear
E_U = m(Uold, Aold, W, vel) \
      - k(Uold, Aold, Phiold, W )

# Weak form for Phi... For the fracture case ----------------------------------#
E_Phi = (-mu2*inner(grad(Phiold),vel)*Q + Gc*l*inner(grad(Phiold),grad(Q)) + ((Gc/l)+2.0*psi_elastic(Uold))*inner(Phiold,Q) - 2.0*psi_elastic(Uold)*Q)*dx

form = E_A + E_U + E_Phi

# Nonlinear Solver
J = derivative(form, old, dU)
problem = NonlinearVariationalProblem(form, old, bc_U, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
info(prm, True)
prm["nonlinear_solver"] = "newton"
prm["newton_solver"]["linear_solver"] = "mumps"
prm["newton_solver"]["absolute_tolerance"] = 1E4
prm["newton_solver"]["relative_tolerance"] = 1E4
prm["newton_solver"]["maximum_iterations"] = 100

#########################################################################################
# Initialization of the iterative procedure and output requests
time = np.linspace(0, T, Nsteps+1)
tol = 1e10 

#Name of file used for storage
store_Phi = File ("mydata/Phi.pvd")
store_u = File ("mydata/u.pvd")
store_a = File ("mydata/a.pvd")
store_stress_elas = File ("mydata/stress_elas.pvd")
store_eta_k = File ("mydata/eta_k.pvd")

#########################################################################################
# Storing things at t = 0:
store_a << Aold
store_u << Uold # mm
store_Phi << Phiold
stress_elas.assign(project( sigma(Uold, Aold, Phiold) ,sigma_fs, solver_type="cg", preconditioner_type="amg")) # 1MPa = 1N/mm^2
store_stress_elas << stress_elas
eta_k.assign( project( eta_k, VV, solver_type="cg", preconditioner_type="amg"))
store_eta_k << eta_k

print ('Saving initial condition')
#########################################################################################
# Looping through time here.
for (i, dt) in enumerate(np.diff(time)):

    t = time[i+1]
    print("Time: ", t)
    iter = 0
    err = 1e11

    while err > tol:
        iter += 1
        # Nonlinear Solver
        solver.solve()
        (Anew, Unew, Phinew) = old.split(True)

        # Calculate error
        err_a = errornorm(Anew,Aold,norm_type = 'l2',mesh = None)
        err_u = errornorm(Unew,Uold,norm_type = 'l2',mesh = None)
        err_Phi = errornorm(Phinew,Phiold,norm_type = 'l2',mesh = None)
        err = max(err_a, err_u, err_Phi)
        print(err_a, err_u, err_Phi)
	
        if err < tol:
            # Update old fields from previous timestep with new quantities
            Aold.assign(Anew)
            Uold.assign(Unew)
            Phiold.assign(Phinew)

            if round(t*1e1) % 2 == 0: # each saved data point is 2e-9s
                store_a << Aold
                store_u << Uold #mm
                store_Phi << Phiold
                stress_elas.assign(project( sigma(Uold, Aold, Phiold) ,sigma_fs, solver_type="cg", preconditioner_type="amg")) # 1MPa = 1N/mm^2
                store_stress_elas << stress_elas

                File('mydata/saved_mesh.xml') << mesh_half
                File('mydata/saved_a.xml') << Aold
                File('mydata/saved_u.xml') << Uold
                print ('Iterations:', iter, ', Total time', t, 'Saving datapoint')
 	    
print ('Simulation completed') 
#########################################################################################

Thank you so much!

Regards
Janel

5

Please start by checking that the boundary conditions are properly applied to a single vector.
Here is an example of how to code that up for your example:

# Janel Chua
#########################################################################################
# Preliminaries and mesh
from dolfin import *
import numpy as np

# Define the corner points and resolution
x0, y0 = -0.5, -1  # Bottom-left corner (start point)
x1, y1 = 4.5, 1.0  # Top-right corner (end point)
nx, ny = 500, 200  # Number of divisions in the x and y directions

# Create the rectangular mesh
mesh_half = RectangleMesh(Point(x0, y0), Point(x1, y1), nx, ny)

#########################################################################################
# Define Function Spaces and mixed elements
FE_A = VectorElement("CG", mesh_half.ufl_cell(), 1)
FE_U = VectorElement("CG", mesh_half.ufl_cell(), 1)  # 1 originally # displacement
FE_Phi = FiniteElement("CG", mesh_half.ufl_cell(), 1)  # phi

Z = FunctionSpace(mesh_half, MixedElement([FE_A, FE_U, FE_Phi]))


#########################################################################################
# Boundary conditions
eps = 1000 * DOLFIN_EPS


def left(x, on_boundary):
    return near(x[0], x0, eps) and on_boundary  # (-0.5, x[1])


def right(x, on_boundary):
    return near(x[0], x1, eps) and on_boundary  # (4.5, x[1])


def imposedPhi(x, on_boundary):
    return near(x[0], x0, eps) and on_boundary


bcright = DirichletBC(
    Z.sub(1).sub(0), Constant(0), right
)  # Right displacement loaded: u_1
bcright2 = DirichletBC(
    Z.sub(1).sub(1), Constant(0), right
)  # Right displacement loaded: u_2
stretch1 = Expression(("0.1"), degree=1)
bcleft = DirichletBC(Z.sub(1).sub(0), stretch1, left)  # Left displacement loaded: u_1
phiright = DirichletBC(Z.sub(2), Constant(0), right)
phileft = DirichletBC(Z.sub(2), Constant(1.0), imposedPhi)

# Dirichlet Boundary Condition used:
bc_U = [bcright, bcright2, bcleft, phiright, phileft]


wh = Function(Z)
wh.vector()[:] = 0.1
[bc.apply(wh.vector()) for bc in bc_U]
print(wh.vector()[:])
for i in range(Z.num_sub_spaces()):
    with XDMFFile(f"wh{i}.xdmf") as file:
        file.write(wh.split()[i])

Now you can visually inspect w0.xdmf, w1.xdmf and w2.xdmf and check if you get the expected boundary condition.

Your code is still far from minimal, and very hard for me to access. You need to go through it step by step, and be much clearer on what works and what doesn’t work.