Boundary conditions on subdomain

Question

Hi, I'd like to specify a boundary condition as being on the edge of a mesh subdomain, rather than having to supply a particular condition in my code (e.g. near(x[0], 1)). I have read the fenics tutorial and have successfully made mesh subdomains (Sec 4.5 in the tutorial vol 1), but am uncertain as to how to specify the edges of these subdomains for boundary conditions.

For example I might want to solve the Laplace equation, del^2 u = 0, for the geometry shown below, with the solution u=5 at the boundaries of the arbitrary shape near the centre and u=0 at the edges of the square domain. How would I specify the edges of the arbitrary shape as a DirichletBC without specifying any coordinates?

Can somebody please post a simple example, and if possible a complete script with simple mesh definitions (e.g. a circle inside a square).

Thanks for your help.

Answer 1

Hi,

I had the same problem. Here is an example that might help you:

from dolfin import *
from mshr import *
from math import pi

# -------------------- #
#     Create  mesh     #
# -------------------- #
rectangle = Rectangle(Point(-1., -1.), Point(1., 1.))
R = 0.5
origin = Point(0.,0.)
circle = Circle(origin, R, segments=32)
domain = rectangle
domain.set_subdomain(1, circle)
domain.set_subdomain(2, rectangle - circle)
mesh = generate_mesh(domain, 15) # 15 is the resolution


# Create subdomains
subdomains = MeshFunction("size_t", mesh, 2, mesh.domains())


# -------------------- #
#   Create boundaries  #
# -------------------- #
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] + 1.) < DOLFIN_EPS

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - 1.) < DOLFIN_EPS

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] + 1.) < DOLFIN_EPS

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - 1.) < DOLFIN_EPS


boundaries = FacetFunction("size_t", mesh, 0 ) # this index 0 is an alternative to the command boundaries.set_all(0)

Bottom().mark(boundaries, 1)
Left().mark(boundaries, 2)
Right().mark(boundaries, 2)
Top().mark(boundaries, 3)

# ---------------------------- #
#   Create boundary interface  #
# ---------------------------- #
# Option 1: 
    for f in facets(mesh):
   domains = []
   for c in cells(f):
       domains.append(subdomains[c])
   domains = list(set(domains)) # remove duplicates
   if len(domains) > 1: # it is on the interface
       boundaries[f] = 4
# Option 2:
#for f in facets(mesh): # For all edges in the mesh
#    p0 = Vertex(mesh, f.entities(0)[0]) # save the two ending points p0 and p1
#    p1 = Vertex(mesh, f.entities(0)[1])
#    if near(p0.x(0)*p0.x(0) + p0.x(1)*p0.x(1), R*R, 0.05) and near(p1.x(0)*p1.x(0) + p1.x(1)*p1.x(1), R*R, 0.05): # check if the points lie on the circle - if yes, put a label on this edge
#        boundaries[f] = 4



# --------------------------------- #
# Chech that everything is correct  #
# --------------------------------- #

# Define a function of ones on the mesh
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.vector()[:] = 1.

# dS on the inteface = 2*pi*R
dS = dS(subdomain_data=boundaries)
form4 = f*dS(4) # integration on the whole facet marked with number 4
print form4
ans = assemble(form4)   
print ans, "vs", 2*pi*R

# dS on an external boundary = 0
form1 = f*dS(1) 
ans = assemble(form1)  
print ans, "vs", 2., "~>", "dS is non zero only in the boundaries' interface!"

# ds on an external boundary = 2
ds = ds(subdomain_data=boundaries)
form1 = f*ds(1) # integration on the external boundary marked with number 1
ans = assemble(form1)   
print ans, "vs", 2.

Comments

Thank you for responding,

If you had the same problem, then how did you resolve the problem? The example you have posted is interesting, and I could see aspects of it as being a part of the solution, but it is also not a direct answer to my question. No boundary conditions are applied in this example. Do you think you could answer the question more directly? How can I use this to apply boundary conditions to these subdomains.

Thanks again.

---

Hi again,

I've given this a bit more of a go and by working from your example I have been successful. This method isn't particularly elegant - it's looping through all facets, as per the suggestion from meron, but it has managed to get the job done. You will notice a part of your example code within the script, thank you.

My final (successful) code is as follows

import fenics as fn
import mshr as ms

# Define Constants
domain_width = 10.0
circ_rad = 2.0

#Define Mesh Objects
domain = ms.Rectangle(fn.Point(-domain_width/2, -domain_width/2),
                    fn.Point(domain_width/2, domain_width/2))
circ = ms.Circle(fn.Point(0.0, 0.0), circ_rad/2)

# Make subdomain from mesh elements
domain.set_subdomain(2, circ)

# Generate Mesh 
mesh = ms.generate_mesh(domain, 200)

# Define boundaries 
boundary_markers = fn.MeshFunction('size_t', mesh, 2, mesh.domains())
boundaries = fn.MeshFunction('size_t', mesh, 1, mesh.domains())

# Use the cell domains associated with each facet to set the boundary
for f in fn.facets(mesh):
    domains = []
    for c in fn.cells(f):
        domains.append(boundary_markers[c])

    domains = list(set(domains))
    if len(domains) > 1:
        boundaries[f] = 2


# Make Function Space
V = fn.FunctionSpace(mesh, 'P', 2)

# Boundary on central shape
box_bc = fn.DirichletBC(V, fn.Constant(5), boundaries, 2)

# Outer Boundary
domain_bc = fn.DirichletBC(V, fn.Constant(0), 'on_boundary')
bcs = [box_bc, domain_bc]

# Solve
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx
L = fn.Constant('0') * v * fn.dx
u = fn.Function(V)
fn.solve(a == L, u, bcs)

# Output
outputFile = fn.File('postOutput.pvd')
outputFile << u

The output of this script is shown in the image below

Note that my boundary is labelled 2, and that this can be changed if the script is to be used by others. I believe that this script could quite easily be modified to, for example, set boundary conditions on two shapes by replacing the

    if len(domains) > 1: # it is on the interface
        boundaries[f] = 2

with something that checks which domains are listed and assigns the appropriate boundary identifier.

I will mark this as the best answer, since I have actually got a working solution from it, but I think it is fair to say that looping through all facets shouldn't be necessary. As meron demonstrated boundaries are not automatically identified, while the subdomains are. This feels like a bug; although I am a new user and it is possible that there's a very good reason for this that I am unaware of. In any case I feel that setting a boundary condition on a geometric object defined with the software should not be so involved.

Thank you very much for your help caterinabig and meron.

---

I agree, it's not very elegant, but it works and you have to do it only once when you generate the mesh.

Indeed, once you know how to mark facets on the interface (in your case with number 2) then you simply use the DirichletBC function as for any other boundary:

box_bc = fn.DirichletBC(V, fn.Constant(5), boundaries, 2)

Answer 2

Hey there,

I think what you need is covered in the tutorial here. Does that help?

Comments

Hi thank you for responding. As far as I can see what I want is not covered in that part of the tutorial. In the linked example the subdomain is specified with

tol = 1E-14
class Omega_0(SubDomain):
        def inside(self, x, on_boundary):
            return x[1] <= 0.5 + tol

however my question is how to avoid having to specify things like x[1] <= 0.5.

Perhaps I could be more clear. If I define a mesh with, for example,

import fenics as fn
import mshr as ms

domain_width = 10.0
circ_rad = 2.0

domain = ms.Rectangle(fn.Point(-domain_width/2, -domain_width/2),
                      fn.Point(domain_width/2, domain_width/2))                                                                      
circ = ms.Circle(fn.Point(0.0, 0.0), circ_rad/2)                                                                                 
domain.set_subdomain(2, circ)
mesh = ms.generate_mesh(domain, 30)

This makes a mesh with a clear circle. I would then like to specify a DirichletBC on the edge of this circle, without having to manually specify the equation of the circle.

I have not yet achieved this, however I have tried things such as

V = fn.FunctionSpace(mesh, 'P', 2) 
boundary_markers = fn.MeshFunction('size_t', mesh, 1, mesh.domains())

box_bc = fn.DirichletBC(V, fn.Constant(5.0), boundary_markers, 2)
domain_bc = fn.DirichletBC(V, fn.Constant(0.0), 'on_boundary')                                                                                                      
bcs = [box_bc, domain_bc]

However when using this to solve, I get "Warning: Found no facets matching domain for boundary condition." and the BC is not applied anywhere. I'd like to know where I am going wrong.

Thanks again for your help.

---

What can go wrong is that with your specified boundary marker, the mesh is too coarse to "see" the circle- one fix is to make your mesh less coarse, which should be straightforward to achieve. Does that work?

---

I have just tried making my mesh more dense by using 200, rather than 30 in the generate_mesh() function. It now reads

mesh = ms.generate_mesh(domain, 200)

This has not fixed the issue, and I still receive the "*** Warning: Found no facets matching domain for boundary condition."

---

Can you find cells in the subdomain defined by the circle?
You can find how to do this here.

If not, the subdomain was not created for some reason.

---

I'm not sure how to do this when I have specified my geometry with mshr, and then haven't manually made a SubDomain child class.

If you know how to achieve what I'm asking for please could you post an example code.

---

Hi again,

I wasn't sure how to check this using the methods on that thread. Their examples had a subdomain being defined as a SubDomain child class. My subdomains are made as a part of the mesh, and thus I don't have a manually defined SubDomain child class.

To try and overcome this, I've tried to tag the cells inside the subdomain. To do this I changed my meshfunction to be 2D and then made an expression, which allows me to check the cell index. I applied a Dirichlet boundary to every point (set with a function that always returns true), with a value set by this expression. This has allowed me to highlight cells which belong to different domains. My full code is

import fenics as fn
import mshr as ms

# Define Constants
domain_width = 10.0
circ_rad = 2.0

#Define Mesh Objects
domain = ms.Rectangle(fn.Point(-domain_width/2, -domain_width/2),
                  fn.Point(domain_width/2, domain_width/2))
circ = ms.Circle(fn.Point(0.0, 0.0), circ_rad/2)

# Make subdomain from mshr geometric elements
domain.set_subdomain(2, circ)

# Generate Mesh and MeshFunction
mesh = ms.generate_mesh(domain, 200)
boundary_markers = fn.MeshFunction('size_t', mesh, 2, mesh.domains())

# Make Function Space
V = fn.FunctionSpace(mesh, 'P', 2)

# Boundary Conditions
class tempClass(fn.Expression):
    '''
    Temporary expression for setting dummy boundary condition
    '''
    def __init__(self, markers, **kwargs):
    self.markers = markers

    def eval_cell(self, values, x, cell):
    if self.markers[cell.index] == 2:
        values[0] = 5.0
    else:
        values[0] = 0.0

def tempFunc(x, on_boundary):
    '''
    Temporary function for setting dummy boundary condition
    '''
    return True

# Boundary on central shape
box_potential = tempClass(boundary_markers, degree=1)  # Expression
box_bc = fn.DirichletBC(V, box_potential, tempFunc)

# Outer Boundary
domain_bc = fn.DirichletBC(V, fn.Constant(0), 'on_boundary')
bcs = [box_bc, domain_bc]

# Solve
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
a = fn.dot(fn.grad(u), fn.grad(v)) * fn.dx
L = fn.Constant('0') * v * fn.dx
u = fn.Function(V)
fn.solve(a == L, u, bcs)

# Output
outputFile = fn.File('postOutput.pvd')
outputFile << u

The output of this is shown in the image below. Clearly the circle is being identified as being in the subdomain of index 2. However I still do not know how to set a Dirichlet boundary condition on the edge of this subdomain in a similar way.

---

Ok, that looks good. So you are able to set boundary conditions in the circle that appears yellow in the plot above, but what you want is to pose the boundary condition only on the outer edge? Clearly, this must have something to do with the way that circ is defined- everything inside the circle has a subdomain marker 2 . A solution is to define a second circle with a radius slightly bigger than the first one and then define Dirichlet BCs on the resulting annulus. E.g. if the second circle has a subdomain marker 4, then use

if self.markers[cell.index] != 2 and self.markers[cell_index] == 4:

to pose a BC on the annulus.

---

Sorry, but the script above is nothing like what I actually want to achieve, it just highlights that the circle is being placed into a subdomain. I have not set a boundary condition on the circle, I "selected" it with an expression, and applied a boundary condition to every single point on the mesh; that's why the solution has every point as either 0 or 5. That's not much better than drawing a circle in mspaint, and it is not the same as applying a boundary condition on the circle.

It makes little difference to me whether I apply the boundary condition just to the edge of the circle or to the whole circle, but I want to apply the BC to the subdomain that the circle is in. Do you know how to do this?

---

Hi!
It seems you can get subdomains (dx measure) but not boundaries (ds/dS) from the mesh.
If you look at the line you used to compute the boundary markers on a facetfunction and inspect them

boundary_markers = fn.MeshFunction('size_t', mesh, 1, mesh.domains())

import numpy as np

print np.unique(boundary_markers.array()[:])

You can indeed see that no boundaries have been set.
When using a cellfunction (change the 1 to 2) you can see the regions, but this can't be used to impose a boundary condition.
Since I don't know how (and whether) the boundaries can actually be found in the mesh, you could write a script to compute the boundaries based on the subdomains (if you want I could send a script). Not a very nice solution, but it would work.