Conditional in an exact solution defined as a UserExpression by subclassing

Dear Colleagues,

I’m applying the exact solution in an infinite plate with a circular inclusion I used the ufl expression conditional (lt(r,R), u_in, u_out) but it show the error division by zero, if someone could help me to know why it shows this error I’ll be grateful , I share my code:

``````domain = Rectangle(Point(-5, -5), Point(5, 5))
mesh = generate_mesh(domain, 50)
V = VectorFunctionSpace(mesh,'P', 1)

r = Expression("hypot(x[0], x[1])", degree = 2)
theta = Expression("atan2(x[1], x[0])", degree = 2)

T = Constant(1.0)
R = Constant(sqrt(2.0)) # radius of the carbon fiber

E_1 = Constant(4.5E9)
E_ = Constant(39E10)
nu_1 = Constant(0.4)
nu_ = Constant(0.35)

lamb_1 = Constant((E_1*nu_1)/((1 + nu_1)*(1 - 2*nu_1)))
mu_1 = Constant(E_1/(2*(1 + nu_1)))
lamb_ = Constant((E_*nu_)/((1 + nu_)*(1 - 2*nu_)))
mu_ = Constant(E_/(2*(1 + nu_)))

k_1 = Constant((3 - nu_1)/(1 + nu_1))
k_ = Constant((3 - nu_)/(1 + nu_))

class Exact_disp(UserExpression):
def __init__(self, T, mu_, mu_1, k_1, k_, R, **kwargs):
super().__init__(**kwargs)
self._T = T
self._mu_ = mu_
self._mu_1 = mu_1
self._k_1 = k_1
self._k_ = k_
self._R = R

def eval(self, value, x):
T = self._T
mu_ = self._mu_
mu_1 = self._mu_1
k_1 = self._k_1
k_ = self._k_
R = self._R
r = hypot(x[0], x[1])
theta = atan2(x[1], x[0])
alpha_ =  2*mu_ - mu_1 + k_*mu_1
beta_ = k_1*mu_ + mu_1
u_r= conditional(lt(r,R), (r*T*(-1 + k_)*(1 + k_1))/(16*mu_ + 8*(k_ - 1)*mu_1), ((T/(4*mu_1))*(pow(R,2)/r + r*0.5*(k_1 - 1))) - (pow(R,2)*T*mu_*(1 + k_1)/(4*r*mu_1*alpha_)) +\
(r*T*cos(2*theta)/(4*mu_1)) + ((pow(R,2)*T*cos(2*theta))/(4*beta_))*(mu_/(mu_1 - 1))*((pow(R,2)/pow(r,3)) - (1 + k_1)/r))
u_t = conditional(lt(r,R), (-0.25*r*T*sin(2*theta)*(1 + k_1))/(k_1*mu_ + mu_1), (T*sin(2*theta)/(4*pow(r,3)*mu_1*beta_))*(pow(R,4)*(-mu_ + mu_1) + pow(r,2)*pow(R,2)*(k_1 - 1)*(mu_1 - mu_) + pow(r,4)*beta_))
value[0] = u_r*cos(theta) - u_t*sin(theta)
value[1] = u_r*sin(theta) + u_t*cos(theta)
def value_shape(self):
return(2,)

u_exact = Exact_disp(T=T, mu_=mu_, mu_1=mu_1, k_1=k_, k_=k_, R=R)

Exact_space = VectorFunctionSpace(mesh, 'P', 1)
u_D = project(u_exact, Exact_space)
``````

the Error:
Traceback (most recent call last):
File “inclusion_2D.py”, line 172, in
u_D = project(u_exact, Exact_space)
File “/usr/lib/python3/dist-packages/dolfin/fem/projection.py”, line 133, in project
form_compiler_parameters=form_compiler_parameters)
File “/usr/lib/python3/dist-packages/dolfin/fem/assembling.py”, line 382, in assemble_system
assembler.assemble(A_tensor, b_tensor)
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 53, in wrapped_eval
self.user_expression.eval(values, x)
File “inclusion_2D.py”, line 94, in eval
(rTcos(2theta)/(4mu_1)) + ((pow(R,2)Tcos(2theta))/(4beta_))(mu_/(mu_1 - 1))((pow(R,2)/pow(r,3)) - (1 + k_1)/r))
ZeroDivisionError: float division by zero

The `pow(r,3)` term in the denominator goes to zero at points close to origin. I changed your domain to points (5,5) and (15,15) and it runs just fine.

this is weird because actually the pow(r,3) is never close to the origin since the condition for that formula is to be greater than r = 1.414214 = sqrt(2) , then i don’t understand why it could be talk about division by zero if we’re never close to zero