How can define the subdomain

mesh
1

I want to create two domains with different elasticity. the main domain is square with length r2=8, the circle radious is r1=3 inside the square. I wrote the following code, but I am not sure this is correct.
%=================================
E1 = 100e3 # Elastic modulus [MPa]
nu1 = 0.3 # Poisson’s ratio
E2 = 200e3 # Elastic modulus [MPa]
nu2 = 0.25 # Poisson’s ratio

dimension of geometry

r1= 3 #inner radius (mm)
a= 8 #outer radius (mm)
#------------------------------------------
x=mesh_coord
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 8.0)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 8.0)
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return near((sqrt(x[0]**2+x[1]**2)<=3),0.0)

class Omega_2(SubDomain):
def inside(self, x, on_boundary):
return near((sqrt(x[0]**2+x[1]**2)>3),8.0)

Initialize sub-domain instances

left = Left()
top = Top()
right = Right()
bottom = Bottom()
Omega_1 = Omega_1()
Omega_2 =Omega_2 ()

constrained displacement boundary

BC_bottom = DirichletBC(V_1, Constant((0.0)), bottom)
BC_left = DirichletBC(V_0, Constant((0.0)), left)
BC = [BC_bottom, BC_left]
#------------------------------------------------------

Initialize mesh function for interior domains

domains = MeshFunction(“size_t”, mesh, mesh.topology().dim())
domains.set_all(0)
Omega_1.mark(domains, 1)
Omega_2.mark(domains, 2)
#----------------------------

Initialize mesh function for boundary domains

boundaries = MeshFunction(‘size_t’, mesh, dim-1)
boundaries.set_all(0)
left.mark(boundaries,0)
bottom.mark(boundaries,1)
right.mark(boundaries,2)
top.mark(boundaries,3)

2

First of all, please use markdown syntax to highlight your code, see: Problem with bool in solve() - #3 by kamensky

Also note that you should make your code complete, i.e. include all definitions so that people can copy the code and reproduce the behavior.

Finally, you should visualize domains in Paraview for verification, by for instance calling File("domains.pvd")<<domains