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 :slight_smile: , 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 :confused: