Trouble marking boundaries of a 2-layered annulus

[avishmj]

I usually create and mark my geometries through GMSH and upload them to my script. However, in order to create an MWE for a problem, I am trying to create a 2-layered annulus geometry in-script. But I ran into an error regarding marking my boundaries.

Here is the script:

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C3 = Circle(Point(0, 0),3)
C_1 = C2 - C1
C_2 = C3 - C2 
C = C_1 + C_2

mesh = generate_mesh(C, 100)

def inner(x, on_boundary):
    return near(x[0]**2 + x[1]**2, 1) and on_boundary
def middle(x, on_boundary):
    return near(x[0]**2 + x[1]**2, 4) and on_boundary
def middle(x, on_boundary):
    return near(x[0]**2 + x[1]**2, 9) and on_boundary

class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
inner().mark(boundary_markers, 1)

This raises the error:

inner().mark(boundary_markers, 1)
TypeError: inner() missing 2 required positional arguments: 'x' and 'on_boundary'

What’s the error in here and how to fix it? What are the arguments that must be passed on to inner()?

[dokken]

@avishmj

def inner(x, on_boundary):
    return near(x[0]**2 + x[1]**2, 1) and on_boundary

These functions should be classes, similar to:

@avishmj

class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True

[avishmj]

Ok. So I changed the code to the following:

class boundary_i(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1) and on_boundary
class boundary_m(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4) and on_boundary
class boundary_o(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 9) and on_boundary

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary_i().mark(boundary_markers, 1)

ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
print(assemble(Constant(1)*ds(1)))

Now the error is

print(assemble(Constant(1)*ds(1)))
  File "/usr/lib/python3/dist-packages/ufl/measure.py", line 429, in __rmul__
    error("This integral is missing an integration domain.")
  File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: This integral is missing an integration domain.

what’s wrong here?

[dokken]

@avishmj

ds = Measure('ds', subdomain_data=boundary_markers)

You need to add:

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

[avishmj]

That works but it’s giving me the wrong results. Here is the code:

C1 = Circle(Point(0, 0),1)
C2 = Circle(Point(0, 0),2)
C3 = Circle(Point(0, 0),3)
C_1 = C2 - C1
C_2 = C3 - C2 
C = C_1 + C_2

mesh = generate_mesh(C, 100)

class boundary_i(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1) and on_boundary
class boundary_m(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 4) and on_boundary
class boundary_o(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 9) and on_boundary

class mesh1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]**2 + x[1]**2 > 1 and x[0]**2 + x[1]**2 < 4
class mesh2(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0]**2 + x[1]**2 > 4 and x[0]**2 + x[1]**2 < 9

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary_i().mark(boundary_markers, 1)
boundary_m().mark(boundary_markers, 2)
boundary_o().mark(boundary_markers, 3)
mesh1().mark(surface_markers, 4)
mesh2().mark(surface_markers, 5)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)
print(assemble(Constant(1)*ds(1)))
print(assemble(Constant(1)*dx(4)))

output:

0.0
0.0

But the command print(assemble(Constant(1)*ds(1))) is supposed to give me the circumference of the circle of radius 1 and print(assemble(Constant(1)*dx(4))) is supposed to give me \pi(2^2-1^2) since dx(4) is the surface measure of the annulus bounded by circles of radii 1 and 2.

What did I do wrong?

[dokken]

I would suggest you output your boundary markers and surface markers to pvd/xdmf and visualize them in paraview, to check if the correct regions are actually marked.

[avishmj]

Neither the boundary markers nor the surface markers are featuring in the pvd output file.

[dokken]

For the surface markers, you should not use the on_boundary argument, as a cell is not on the boundary. I would suggest the following:

class mesh2(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 > 4
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 4)
mesh2().mark(surface_markers, 5)

For the boundary markers I would use something like:

class boundary_i(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0]**2 + x[1]**2, 1, eps=1e-3) and on_boundary


class boundary_o(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0]**2 + x[1]**2, 9, eps=1e-3) and on_boundary


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
boundary_i().mark(boundary_markers, 1)
boundary_o().mark(boundary_markers, 3)

I am not sure how you want the interior marking (boundary_m) to work, as this interface is not aligning with the mesh.

[avishmj]

@dokken

I am not sure how you want the interior marking (boundary_m) to work, as this interface is not aligning with the mesh.

I do need the middle boundary in order to define two different domains, the outer annulus and the inner one to define two different functions on the two layers.
I am actually trying to recreate an error that I am experiencing with a three-layer annulus. Here is the problem/error I am trying to reproduce:

For a three layer geometry created in GMSH, like

when I am creating three functions psi1, psi2, and psi3 on the surfaces tagged 5,6, and 7 and solving a variational problem with the function

Pi = psi1*dx_int(5) + psi2*dx_med(6) + psi3*dx_ad(7) - dot(-T*(n), u)*ds_int(1)

through the commands

F = derivative(Pi, u, v)

# Compute Jacobian of F
Ju = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcs, J=Ju,
form_compiler_parameters=ffc_options)

FEniCS is returning the error:

solve(F == 0, u, bcs, J=Ju,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 308, in _solve_varproblem
    problem = NonlinearVariationalProblem(eq.lhs, u, bcs, J,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/problem.py", line 157, in __init__
    F = Form(F, form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 23, in __init__
    self.subdomains, = list(sd.values())  # Assuming single domain
ValueError: too many values to unpack (expected 1)

Now I don’t understand the source of the error. So I am trying to reproduce the error with a two-layer annulus or want to check if the error still persists when I am defining the mesh in-script.

[avishmj]

Ok. I got it. Problem solved. The problem was that while defining Pi, I used the submesh measures and not the global measure.

[dokken]

Btw, see Solving PDEs in Python - <br> The FEniCS Tutorial Volume I on how to mark subdomains with mshr

[avishmj]

@dokken

Btw, see Solving PDEs in Python -
The FEniCS Tutorial Volume I on how to mark subdomains with mshr

I have a couple of questions. One regarding one of the features I saw in the link you suggested and also in other examples, and one regarding the solution to the problem I posted.

1. I should always define my subdomains before I generate the mesh for the whole domain. In that way, the boundaries of the submesh will align with the mesh nodes (i.e., the subdomain boundary will not intersect a cell in a way that none of the nodes of the cell lies on the boundary). On the other hand, if I generate a mesh and then define my subdomain, it may happen that the subdomain boundaries intersect a cell with none of the nodes of that cell being on the subdomain boundary. Is that correct? And is that what you meant when in a previous post you mentioned that the subdomain boundary (boundary_m) doesn’t align with the mesh?

2. In the solution to the problem I posted before your last comment, why are the subdomain measures not working when I am defining the function Pi? Why do I need to plug in the global measures that I defined when I imported the mesh and before extracting the subdomains? According to me, subdomain measures are subsets of the global measures. So they should both work equally well. Here is an explanation that one of my friends suggested and I am confused as to why is that true.
   “The reason for using submeshes is to solve a problem independent of the others. However, if you want to integrate over the whole domain with different subdomain integrations you just label the global measure.”

[dokken]

@avishmj

I should always define my subdomains before I generate the mesh for the whole domain. In that way, the boundaries of the submesh will align with the mesh nodes (i.e., the subdomain boundary will not intersect a cell in a way that none of the nodes of the cell lies on the boundary). On the other hand, if I generate a mesh and then define my subdomain, it may happen that the subdomain boundaries intersect a cell with none of the nodes of that cell being on the subdomain boundary. Is that correct? And is that what you meant when in a previous post you mentioned that the subdomain boundary (boundary_m) doesn’t align with the mesh?

Yes, that is what I meant.

@avishmj

In the solution to the problem I posted before your last comment, why are the subdomain measures not working when I am defining the function Pi? Why do I need to plug in the global measures that I defined when I imported the mesh and before extracting the subdomains? According to me, subdomain measures are subsets of the global measures. So they should both work equally well. Here is an explanation that one of my friends suggested and I am confused as to why is that true.
“ The reason for using submeshes is to solve a problem independent of the others. However, if you want to integrate over the whole domain with different subdomain integrations you just label the global measure. ”

I cannot answer this question, as you in the previous question only added one line of your code, and not a minimal example. Thus I do not know how you are defining your measures.