Boundary Conditions not working for mixed element method on a nonlinear problem

General
1

Greetings,

I am running into problems applying boundary conditions when using the mixed element method in FEniCS. Specifically, there are no runtime errors but when i check the results, the boundary conditions like phicrack and bcleft1, bcleft21 are simply not applied.

Here is the segment on how i am applying dirichlet boundary conditions:

# Boundary conditions
def left(x,on_boundary):
    return near(x[0],-0.5) and on_boundary #(-0.5, x[1])
def leftTopHalf(x, on_boundary):
    return near(x[0],-0.5) and (x[1] > 0) and on_boundary #(-0.5, x[1])
def leftBotHalf(x, on_boundary):
    return near(x[0],-0.5) and (x[1] < 0) and on_boundary #(-0.5, x[1])
def right(x,on_boundary):
    return near(x[0],5) and on_boundary #(5, x[1])
def imposedPhi(x,on_boundary):
    return abs(x[1]) < 2e-02 and x[0] <= 0 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.01"),degree=1)
stretch2 = Expression(("-0.01"),degree=1)
bcleft1= DirichletBC(Z.sub(1).sub(1), stretch1, leftTopHalf) #Left displacement loaded: u_1
bcleft21 = DirichletBC(Z.sub(1).sub(1), stretch2, leftBotHalf) #Left displacement loaded: u_2

phiright = DirichletBC(Z.sub(2), Constant(0), right)
phicrack = DirichletBC(Z.sub(2), Constant(1.0), imposedPhi)

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

And below is the minimal working code (i tried to minimize it as much as possible sorry if its still too long)

# 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')

# 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

#########################################################################################
# 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) and on_boundary #(-0.5, x[1])
def leftTopHalf(x, on_boundary):
    return near(x[0],-0.5) and (x[1] > 0) and on_boundary #(-0.5, x[1])
def leftBotHalf(x, on_boundary):
    return near(x[0],-0.5) and (x[1] < 0) and on_boundary #(-0.5, x[1])
def right(x,on_boundary):
    return near(x[0],5) and on_boundary #(5, x[1])
def imposedPhi(x,on_boundary):
    return abs(x[1]) < 2e-02 and x[0] <= 0 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.01"),degree=1)
stretch2 = Expression(("-0.01"),degree=1)
bcleft1= DirichletBC(Z.sub(1).sub(1), stretch1, leftTopHalf) #Left displacement loaded: u_1
bcleft21 = DirichletBC(Z.sub(1).sub(1), stretch2, leftBotHalf) #Left displacement loaded: u_2

phiright = DirichletBC(Z.sub(2), Constant(0), right)
phicrack = DirichletBC(Z.sub(2), Constant(1.0), imposedPhi)

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

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') 
#########################################################################################

Trial14_coarseEdges8.xml is just a 2D rectangular domain with a physical notch at y=0 going from -0.5<x<0 (where the phi boundary condition is applied). The Domain ranges from -0.5<x<5 and -1<y<1

Really appreciate any suggestions! Thank you so much!

2

Are none of the boundary conditions being imposed? What have you done to debug this issue? Have you visualized your boundary equations (to make sure they are actually being applied) by writing a FacetFunction to file? Without the mesh it’s difficult to help but you use near() which without a third argument defaults to DOLFIN_EPS.

1 Like
3

Hi Sven,

Thank you so much for replying me! none of the dirichlet BC are being imposed. I checked the Phi outputs on paraview to check if phileft was imposed (if it is there should be a red line on the left boundary) and it wasn’t. Neither was bcleft imposed (i displaced the left boundary by 0.1 but when i turn on warp by vector there was no displacement at all).

Here is an edited code with the mesh :slight_smile:

# 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
# 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

#########################################################################################
# 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) and on_boundary #(-0.5, x[1])
def right(x,on_boundary):
    return near(x[0],4.5) and on_boundary #(4.5, x[1])
def imposedPhi(x,on_boundary):
    return near(x[0],-0.5) 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') 
#########################################################################################

Appreciate your help very much!

Regards
Janel