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