Dirichlet BC's and mesh subdomain in 2D

Question

Dear all,

I wonder if there is a specific way for defining Dirichlet BC's when working with mesh subdomains. I use the following script to generate a 2D mesh with mshr (file meshgen.py):

import dolfin
from mshr import *
def set_mesh(ref,h,w):
   p1 = dolfin.Point(0.,0.)
   p2 = dolfin.Point(w,h)
   p3 = dolfin.Point(w,1/2.0*h)
   domain =   Rectangle(p1,p2)
   domain.set_subdomain(1, Rectangle(p1, p3))
   mesh2d = generate_mesh(domain, ref)
   return (mesh2d)

And I noticed some difficulties when I apply Dirichlet BC's with the "standard" procedure :

from meshgen import set_mesh
from dolfin import *

mesh= set_mesh(10,1.,1.2)

V = VectorFunctionSpace(mesh, "Lagrange", 1)

class Corner(SubDomain):
  def inside(self, x, on_boundary):
    return x[0]==0.0 and x[1]==0.0

class Border(SubDomain):
  def inside(self, x, on_boundary):
    return x[1]==0.0

border = Border()
corner = Corner()
bc1 = DirichletBC(V, Constant((0.,0.)),corner,method="pointwise")
bc2 = DirichletBC(V.sub(1), 0., border)
bcs = [bc1, bc2]

Is there a specific way to define BC's when working with mesh subdomains ?

thank you in advance,

Answer 1

Hi, it is not clear what your problem is but I suspect you get warnings about bc2 being unable to find facets to apply boundary conditions at. This is caused by the definitions of your subdomains where you use == to compare two floats. You should use instead the near function which checks equality of its arguments within some tolerance. That is, your domains could look as follows

TOL = DOLFIN_EPS  # Or something larger 

class Corner(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[0], 0., TOL) and near(x[1], 0., TOL)

class Border(SubDomain):
  def inside(self, x, on_boundary):
    return near(x[1], 0., TOL) 

Comments

Yes correct, and it works so !

I also found an alternative solution using gmsh to create the mesh, and then using meshconvert:

from dolphin_utils.meshconvert import meshconvert
meshconvert.convert2xml(subdir + meshname + ".msh", subdir + meshname + ".xml","gmsh")

I can assign indicators to boundary parts directly in the .geo file, for instance :

Physical Line(%) = {%};

Thanks for your answer !