Interior boundary conditions

Question

Dear all,

I am trying to solve a boundary problem. As I understand, one way to do that is to apply "boundary conditions" on the interior of the domain, somehow in the opposite way as boundary conditions are commonly used. I first try on vectors before solving a system. I face the problem that the boundary condition that I declare on the interior of the domain (i.e. "not on_boundary") are also applied on the domain boundary. Am I doing something wrong ?

Here is my test code:

from dolfin import *

mesh = RectangleMesh(0, 0, 1, 1, 10, 10)

V = FunctionSpace(mesh, "CG", 1)

def inside(x, on_boundary):
    return not on_boundary

u_in = Constant(-1.0)

bc_in  = DirichletBC(V, u_in,  inside)

u = Function(V)
u.vector()[:] = 0.5
bc_in.apply(u.vector())

plot(u, axes=True)
interactive()

Here is the output:

I would have expected the domain boundary to be 0.5 everywhere, but apart from 2 corners it has taken the -1.0 value of the interior...

Answer 1

It looks like it's marking all interior facets for application of the DirichletBC. All nodes except for two corner nodes are attached to some interior facet, and hence have the bc applied.

Comments

@tlecomte: Did you try a behaviour of DirichletBC(V, u_in, inside, method='geometric')?

---

I have just tried (I am sorry that I did not see your message ealier), and with DirichletBC(V, u_in, inside, method='geometric') I get the very same result.

---

There is probably a problem that algorithms inside DirichletBC and/or SubDomain::inside does not ensure that not on_boundary is just complement to on_boundary. You can file an issue on bitbucket.

But I would suggest to use rather DirichletBC constructor based on FacetFunction. You can do something like

facet_domains = FacetFunction('size_t', mesh)
#facet_domains.set_all(0) # needed in older versions of DOLFIN
for f in facets(mesh):
    if any(ff.exterior() for ff in facets(f)):
        facet_domains[f] = 666

bc_in = DirichletBC(V, u_in, facet_domains, 0)

---

You could also do something like:

def inside(x, on_boundary):
  return x[0] > DOLFIN_EPS and x[0] - 1.0 < DOLFIN_EPS \
     and x[1] > DOLFIN_EPS and x[1] - 1.0 < DOLFIN_EPS
bc_in = DirichletBC(V, u_in, inside, method = "pointwise")

although this is obviously impractical for more complex geometries.

---

@James: Thanks, but I precisely need to handle complex geometries such as produced with PolygonalMeshGenerator, so that I cannot define a domain using coordinates.

---

I precisely need to handle complex geometries

Then use FacetFunction framework I suggested above. Note that this is very customizable as you can provide your own criterium in a place of any(ff.exterior() for ff in facets(f)).

---

@Jan: Yes, it works beautifully, thank you very much ! here is the output:

---

For completeness, here is the working script (also notice the call to mesh.init() before calling the .exterior() methods):

from dolfin import *

mesh = RectangleMesh(0, 0, 1, 1, 10, 10)

V = FunctionSpace(mesh, "CG", 1)

# initialize mesh connectivity so that Facet.exterior() works
mesh.init()

# define the interior of the domain by looking at each facets
facet_domains = FacetFunction('size_t', mesh)
facet_domains.set_all(0)
for f in facets(mesh):
    if any(ff.exterior() for ff in facets(f)):
        facet_domains[f] = 1

u_in = Constant(-1.0)

bc_in = DirichletBC(V, u_in, facet_domains, 0)

u = Function(V)
u.vector()[:] = 0.5

bc_in.apply(u.vector())

plot(u, axes=True)
interactive()

---

Hello tlecomte,

Can you pls. shed some light on how the below code behaves, especially the for loop.

facet_domains = FacetFunction('size_t', mesh)
facet_domains.set_all(0)
for f in facets(mesh):
 if any(ff.exterior() for ff in facets(f)):
      facet_domains[f] = 1