Conditional in an exact solution defined as a UserExpression by subclassing

1

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

2

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.

3

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: