Define subdomain using other subdomains

mesh
1

Hello,
I am wondering how one might define a subdomain in terms of other subdomains. In the MWE below, the boundary marked as “Rest” does not contain facets directly adjacent to the Left wall.

from fenics import *

mesh = UnitSquareMesh(3,3)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and on_boundary
    
class Rest(SubDomain):
    def inside(self, x, on_boundary):
        return not Left.inside(self,x,on_boundary) and on_boundary
    
domains     = MeshFunction('size_t', mesh, mesh.topology().dim(), 0) 
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)

Left().mark(domains,1) ; Left().mark(boundaries,1)
Rest().mark(domains,2) ; Rest().mark(boundaries,2)

File('boundaries.pvd') << boundaries # save PVD for visualization
2

I’m not exactly sure what your question is here. In your MWE you propose finding “Rest” by taking everything not in “Left” as members of “Rest”. This works just fine, the boundaries.pvd file shows this is correct. Your domains will not have any markers because of your use of “on_boundary” in the class definition of “Rest” and “Left”.

If you are asking for something more complicated where you have many subdomains that are too complicated to be expressed as a class definition (as you’ve provided above), please provide a more specific example. Or, are you asking for other ways to accomplish what your code already does?

3

What I am trying to achieve is for “Rest” to mark all boundaries that are not “Left”. The code in the initial post does not correctly account for facets that are connected to the Left boundary. See the figure below. Notice that the top and bottom left are labelled 0 instead of 2 as desired.

subdomain_mark
subdomain_mark826×697 5.06 KB

Futhermore, we can test that subdomains are marked correctly by integrating over the length of each.

def Edge_Length(ID):
    V = FunctionSpace(mesh, "CG", 1)
    v = Function(V)
    v.assign(Constant(1))
    length = assemble(v*ds(ID))
    return length

print('Left =',Edge_Length(1))
print('Rest =',Edge_Length(2))

which returns,

Left = 1.0
Rest = 2.333333333333333

confirming that the two line segments (each of length 1/3 for this mesh) are missing from ds(2).

My guess is that this is inherent behaviour of the .inside() method: If the segment contains a point on Left it cannot be part of Rest.

So my question is, is there a way to define a subdomain in terms of other subdomains like this?

4

Thank you for the clarification, that was a misunderstanding on my part. I think your assessment is right given your approach “like this”.

I’m not quite sure how to implement this just yet but my initial thought is to get the connected vertices to a given vertex and see if those vertices share multiple exterior facets. Something like this…

class Rest(SubDomain):
    def inside(self, x, on_boundary):
        out = True
        if (Left.inside(self,x,on_boundary)): # If it does lie on another SubDomain
           # Get connected vertices to x
               # for vertex in vertices on_boundary
                   # if all connected vertices are on Left()
                     # out = False

        return out and on_boundary

There might be another way but at that point I think it’s easier in almost all scenarios to just define a second SubDomain explicitly. Especially as the number of SubDomains grows.

5

I would suggest using a slightly different approach.
If you want to have a marker for all entities on the boundaries, start by marking them with a given value (say 2), and then use the other markers to overwrite for the other facets. Minimal example is supplied below:

from fenics import *

mesh = UnitSquareMesh(3,3)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and on_boundary
    
class Rest(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
    
domains     = MeshFunction('size_t', mesh, mesh.topology().dim(), 0) 
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
Rest().mark(domains,2) ; Rest().mark(boundaries,2)
Left().mark(domains,1) ; Left().mark(boundaries,1)


File('boundaries.pvd') << boundaries # save PVD for visualization

I would probably rename Rest to All if I were you.

1 Like
6

Clever. That allows me to define the measures for ds correctly, but what if I wanted to apply a Dirichlet condition on the Rest boundary? As in, DirichletBC(Q, Constant(0), Rest()).

7

How about

DirichletBC(Q, Constant(0), boundaries, 2)
8

That will work. Thanks @dokken!