Dirichlet BC on interior (subdomain) boundary

On Thu, 02 May 2013 14:11:31 -0000
Hossein Nili <email address hidden> wrote:
> New question #228021 on FEniCS Project:
> https://answers.launchpad.net/fenics/+question/228021
>
> Hi
>
> I would much appreciate if you could help me with this:
>
> I am defining subdomains within my main problem domain, and wish to
> impose boundary conditions on the top domain boundary (which I can do
> using straightforward syntax) and also on an interior boundary that
> separates two of the subdomains. When I do the latter using the same
> syntax as used for the former I get this error message:
>
> Warning: Found no facets matching domain for boundary condition.
>
> Is there a different syntax/method for defining Dirichlet BCs on
> interior boundaries?

No. It should matter if boundary is interior or exterior. Can you post
very short complete code?

Jan

>
> Many thanks,
>

Many thanks for responding. I have pasted the code part below. I actually tried doing the same with a very simple script from the FEniCS tutorial, and got the same error message. Maybe I shouldn't be using "on_boundary" in my defining of the boundary conditions on interior boundaries (?)

h = 0.1; hb = 0.1
htopr = h + hb
htopv = 1.01*(h + hb)
mesh = Rectangle(0, 0, 1, htopv, 101, 101)

# Define a MeshFunction over two subdomains
subdomains = MeshFunction('uint', mesh, 2)

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] <= h else False

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] >= h and x[1]<= htopr else False

class Omega2(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[1] >= htopr else False
# note: it is essential to use <= and >= in the comparisons

# Mark subdomains with numbers 0 and 1
subdomain0 = Omega0()
subdomain0.mark(subdomains, 0)
subdomain1 = Omega1()
subdomain1.mark(subdomains, 1)
subdomain2 = Omega2()
subdomain2.mark(subdomains, 2)

V0 = FunctionSpace(mesh, 'DG', 0)
m = Function(V0)

# Loop over all cell numbers, find corresponding
# subdomain number and fill cell value in m
m_values = [1, 0, 0] # values of m in the two subdomains
for cell_no in range(len(subdomains.array())):
    subdomain_no = subdomains.array()[cell_no]
    m.vector()[cell_no] = m_values[subdomain_no]

# Much more efficient vectorized code
# (subdomains.array() has elements of type uint32, which
# must be transformed to plain int for numpy.choose to work)
help = numpy.asarray(subdomains.array(), dtype=numpy.int32)
m.vector()[:] = numpy.choose(help, m_values)

VS = FunctionSpace(mesh, 'Lagrange', 1)

tol = 1E-14

class TopBoundaryV(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - htopv) < tol

class TopBoundaryR(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[1] - htopr) < tol

bc1S = DirichletBC(VS, Constant(1.0), TopBoundaryV())
bc2S = DirichletBC(VS, Constant(1.0), TopBoundaryR())

bcS = [bc1S, bc2S]

On Thu, 02 May 2013 14:36:13 -0000
Hossein Nili <email address hidden> wrote:
> Question #228021 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/228021
>
> Status: Answered => Open
>
> Hossein Nili is still having a problem:
> Many thanks for responding. I have pasted the code part below. I
> actually tried doing the same with a very simple script from the
> FEniCS tutorial, and got the same error message. Maybe I shouldn't be
> using "on_boundary" in my defining of the boundary conditions on
> interior boundaries (?)

You're right. Please, try it first without on_boundary.

Generally, more robust solution is if you're able to pick facets to be
applied bcs on by FacetFunction.

>
>
> h = 0.1; hb = 0.1
> htopr = h + hb
> htopv = 1.01*(h + hb)
> mesh = Rectangle(0, 0, 1, htopv, 101, 101)
>
> # Define a MeshFunction over two subdomains
> subdomains = MeshFunction('uint', mesh, 2)
>
> class Omega0(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] <= h else False
>
> class Omega1(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] >= h and x[1]<= htopr else False
>
> class Omega2(SubDomain):
> def inside(self, x, on_boundary):
> return True if x[1] >= htopr else False
> # note: it is essential to use <= and >= in the comparisons
>
> # Mark subdomains with numbers 0 and 1
> subdomain0 = Omega0()
> subdomain0.mark(subdomains, 0)
> subdomain1 = Omega1()
> subdomain1.mark(subdomains, 1)
> subdomain2 = Omega2()
> subdomain2.mark(subdomains, 2)
>
> V0 = FunctionSpace(mesh, 'DG', 0)
> m = Function(V0)
>
>
> # Loop over all cell numbers, find corresponding
> # subdomain number and fill cell value in m
> m_values = [1, 0, 0] # values of m in the two subdomains
> for cell_no in range(len(subdomains.array())):
> subdomain_no = subdomains.array()[cell_no]
> m.vector()[cell_no] = m_values[subdomain_no]
>
> # Much more efficient vectorized code
> # (subdomains.array() has elements of type uint32, which
> # must be transformed to plain int for numpy.choose to work)
> help = numpy.asarray(subdomains.array(), dtype=numpy.int32)
> m.vector()[:] = numpy.choose(help, m_values)
>
>
> VS = FunctionSpace(mesh, 'Lagrange', 1)
>
> tol = 1E-14
>
> class TopBoundaryV(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[1] - htopv) < tol
>
> class TopBoundaryR(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[1] - htopr) < tol
>
> bc1S = DirichletBC(VS, Constant(1.0), TopBoundaryV())
> bc2S = DirichletBC(VS, Constant(1.0), TopBoundaryR())
>
> bcS = [bc1S, bc2S]
>

I deleted "on_boundary" and it works!

Thank you,

Hossein

An interesting observation is that you'll get the right solution from your equation *only* if the inner boundary on which you are imposing a BC falls along a mesh line. I get substantially different outputs if this rule is not followed: and with irregular inner boundary shapes this is impossible to do.

I am switching to an irregularly-shaped inner boundary, and I think I should doubt the results I get then.

So it looks to me as if "on_boundary" meaning *all over the mesh* is in fact assumed inside the solver, even if (like I did) you omit the "on_boundary" statement: the error message disappears but something unwanted is happening inside disrupting your solution - there should be a way round this, right?

Thanks. I think the pointwise method should handle the problem I am having with irregular boundaries (I was referring to the difficulty in having interior boundaries 'align' with mesh lines - with "pointwise" I believe I can have them follow points rather than lines).

I have a second (related) question which I think is best if I post as a separate question on the forum.

Many thanks for your responses, much appreciated.

Thanks Jan Blechta, that solved my question.