### Defending sub domains

249
views
0
5 months ago by
Dear all,

I am new FEniCS and I am trying to solve simple example of Possion equation. I tried to define sub-domains to apply boundary conditions. I followed the steps in the manual. I have 2 ways to do this so far, the first one is through the setting sub-domains through the main domain and then defining a function to apply the boundary conditions, here is my first code:

domain2 = Rectangle(Point(0, 0), Point(5, 4))

domain2.set_subdomain(1, (Rectangle(Point(1, 2), Point(2, 3)) - Circle(Point(1, 2), 0.9)))
domain2.set_subdomain(2, (Rectangle(Point(3, 2), Point(4, 3)) - Circle(Point(4, 2), 0.9)))
domain2.set_subdomain(3, (Rectangle(Point(1, 0.5), Point(4, 0.75))))

mesh2 = generate_mesh(domain2, 100)

def lowerplate(x):
return (x[0] > 1 + DOLFIN_EPS and x[0] < 4 - DOLFIN_EPS
and x[1] > 0.5 + DOLFIN_EPS and x[1] < 0.75 - DOLFIN_EPS)

def upperplate(x):
return (x[0] > 1 + DOLFIN_EPS and x[0] < 2 - DOLFIN_EPS
and x[1] > 2 + DOLFIN_EPS and x[1] < 3 - DOLFIN_EPS
and (x[0] - 1)**2 + (x[1] - 2)**2 >= (0.9)**2 + DOLFIN_EPS
or (x[0] > 3 + DOLFIN_EPS
and x[0] < 4 - DOLFIN_EPS
and x[1] > 2 + DOLFIN_EPS and x[1] < 3 - DOLFIN_EPS
and (x[0] - 4)**2 + (x[1] - 2)**2 >= (0.9)**2 + DOLFIN_EPS))

V = FunctionSpace(mesh2, "Lagrange", 2)
plot(mesh2)
plt.show()​

The 2nd way to do this is to use the MeshFunction as the manual does, here is my code for it:
domain2 = Rectangle(Point(0, 0), Point(5, 4))

mesh2 = generate_mesh(domain2, 100)

sub_domains = MeshFunction("size_t", mesh2, mesh2.topology().dim())
sub_domains.set_all(0)

class lower_plate(SubDomain):
def inside(self, x, on_boundary):
return True if (x[0] > 1 + DOLFIN_EPS and x[0] < 4 - DOLFIN_EPS
and x[1] > 0.5 + DOLFIN_EPS and x[1] < 0.75 - DOLFIN_EPS) else False

class upper_plate(SubDomain):
def inside(self, x, on_boundary):
return True if (x[0] >= 1 + DOLFIN_EPS and x[0] <= 2 - DOLFIN_EPS
and x[1] >= 2 + DOLFIN_EPS and x[1] <= 3 - DOLFIN_EPS
and (x[0] - 1)**2 + (x[1] - 2)**2 >= ((0.90) + DOLFIN_EPS)**2
or (x[0] >= 3 + DOLFIN_EPS and x[0] <= 4 - DOLFIN_EPS
and x[1] >= 2 + DOLFIN_EPS and x[1] <= 3 - DOLFIN_EPS
and (x[0] - 4)**2 + (x[1] - 2)**2 >= ((0.90) + DOLFIN_EPS)**2)) else False

subdomain_1.mark(sub_domains, 1)
subdomain_2.mark(sub_domains, 2)
V = FunctionSpace(mesh2, "Lagrange", 2)
​

For the first method, I get a really resolved mesh, while this is not the case for the 2nd method (even if I add mesh refinement through MeshFunction, the features are not resolved at all). Here is what I get for the first method (Part of the full domain):

And here is what I get for the 2nd method:

I know that the 2nd figure does not represent the mesh, but even if I added the mesh refinement, The first method is much better (or at least that is what I think)

Can someone explain what is happeneing here ( maybe it is logical error in my code?) Is there a better way to do this? Also, if I want to apply this to 3D will the process be the same?

Thank you for your help
Community: FEniCS Project

2
5 months ago by
Hi,

In your first method, you first defined the subdomain and its boundaries and then generated the mesh. So the mesh generator knows your subdomains.

However, in your second method, the mesh had already been defined and then youv assigned each cell to a subdomain. In this second case, you don't have a mesh that conforms to your subdomains.

It's the same in three dimensions.
1
Oh, I was trying to make a 3d subdomain and got the following message:

*** Error: Unable to setting subdomain.
*** Reason: Subdomains are currently supported only in 2D.
*** Where: This error was encountered inside CSGGeometry.cpp.
*** Process: 0
*** DOLFIN version: 2018.1.0.dev0

3D subdomains are not supported yet. An alternative would be to model the geometry and generate the mesh outside fenics, and then import the mesh with marked subdomains. For example, using gmsh.
written 5 months ago by Igor Baratta