Conditional between a UFL and number

General
1

In a plane strain linear elastic model to calculate the safty factor (SF) in underground mining, i need to implementate a conditional between a UFL and a number to evite obtain NAN results, but it seems like im doing something wrong.

The SF depends from te principal stress values and rock mass propierties, the code is:

#Malla
from dolfin import *
mesh = Mesh("fantasma1.xml")

#Inputs

E = Constant(46.9e3)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
x_max = 1669.4431
y_r=1078.926

#Functions

def eps(v):
  return sym(grad(v))
def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

#Neumann BC
rigth_function = Expression(("-0.027*1.5*(y_r-x[1])", 0), degree=1,y_r=y_r) #Compression in rigth Boundary
tol=10e-10
rigth = AutoSubDomain(lambda x: near(x[0],x_max,tol))
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
rigth.mark(boundaries, 2)
ds = ds(subdomain_data=boundaries)

rho_g = 0.027
f = Constant((0, -rho_g))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)

#Dirichlet BC

def bottom(x, on_boundary):
  return near(x[1], 0.,tol)
def left(x, on_boundary):
  return near(x[0], 0.,tol)
def right(x, on_boundary):
  return near(x[0], x_max,tol)  
bc3 = DirichletBC(V.sub(0), Constant(0.), right)
bc2 = DirichletBC(V.sub(0), Constant(0.), left)
bc1 = DirichletBC(V.sub(1), Constant(0.), bottom)
bc=[bc1,bc2]

#Weak formulation

a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx +inner(rigth_function, u_)*ds(2)
u = Function(V, name="Displacement")
solve(a == l, u, bc)

#Principal Stress

Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u), Vsig)); 

sx = sig.sub(0)
sy = sig.sub(3)
txy = sig.sub(1)

def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2));

S1 = sigma1(sx,sy,txy);
S3 = sigma3(sx,sy,txy)

Vp = TensorFunctionSpace(mesh, "DG",degree=1)
sp = Function(Vsig, name="Principal_Stress")
sp.assign(project(as_tensor(((S1,0.),(0.,S3))), Vp));

#FS

mb=5.157
s=0.02
a=0.502
Em=41.5
sci=65
RCU=296
s3m=sp.sub(3)
s1m=sp.sub(0)

def hb(s3m):
  return(s3m+RCU*(mb*(s3m/RCU)+s)**a)
s1hb= hb(s3m)

a1=abs(s3m)
a2=abs(sci*s**a)

FS= conditional(gt(a1,a2),sci*s**a/s3m,s1hb/s1m)
Vs = TensorFunctionSpace(mesh,"DG",degree=1)
s = Function(Vs,name="FS")
s.assign(project(as_tensor(((FS,0.),(0.,0.))), Vs))

file_results = XDMFFile("fantasma1_Prueba3.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)
file_results.write(sp,0.)
file_results.write(s,0.)

I really hope you can help me

2

Without a minimal reproducible example (on say a built in cube/square) it is really hard to give you any guidance as to what goes wrong, as no-one can run your code.