defining a boundary, w.r.t. a field value

Question

Hello there,
I am quite new to FEniCS, and python too, and I have a question.
Is it possible to define a boundary condition, not only with a geometrical aspect (i.e. x[0] > 0.1), but also with respect to a field value on the boundary (a_field > a_value)?

I gave a naive try, and of course it does not work, saying that there is a type problem when calling the boundary function. And I am afraid I do not know what to do ... Is there something I am missing ?

here is the code :

 from dolfin import *

class boundary(SubDomain):
  def __init__(self,h0):
    self.h0=h0
  def inside(self,x,on_boundary):   
    return (x[0] < 0.1)and(self.h0.vector()>0.1) and on_boundary

#mesh initialisation
km = 1000.
L =160.*km          # length of the domain

N=100               #number of cell in 1 direction
mesh = UnitSquareMesh(N,N)
mesh.coordinates()[:]=mesh.coordinates()*L
deltax = L/N            


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

# setting the topology
Bed= interpolate(Expression("3000+3000*(-exp(-pow(((x[0]-0.6*l)/1.5*0.00003),6))*exp(-pow(((x[1]-0.5*l)/1.5*0.00003),4)) -exp(-pow(x[1]-0.5*l,2)/1.5*0.00000005)*(0.35*exp(-pow(x[0]/1.5*pow(10,-4.8),4))+0.0*exp(-pow(x[0]-0.35*l,2)/1.5*0.00000001)+0.1*exp(-pow(x[0]-0.25*l,2)/1.5*0.00000001))) ",l=L),V)
h0= project(100.,V)
h0.vector()[:]= h0.vector().array() * (Bed.vector().array() <2000.)
Dh0= interpolate(Expression("500*exp(-pow(x[1]-0.5*l,2)/1.5*0.00000005)*(exp(-pow((x[0]+0.1*l)/1.5*pow(10,-4.3),4)))",l=L),V)
Dh0.vector()[:] =Dh0.vector().array() * (Bed.vector().array() < 2500.)

# only a small part reaches the boundary with non zero values (watch the plots) 
h0.vector()[:]=h0.vector() +Dh0.vector()

plot(Bed,title='Bed')
plot(h0,title='h0')
plot(h0+Bed, title='surface')

# should create the Dirichlet BC, but does not ...
BC=boundary(h0)
bc= DirichletBC(V,max(h0.vector()),BC)

#this function is made just to apply the dirichlet BC
bla=interpolate(Expression("0.000001*x[0]"),V)
bc.apply(bla.vector())

interactive()

I would be really gratefull if someone can give me a hint

Comments

You need to be careful that you call the superconstructor when initialising your boundary object. Also, performing the comparison self.h0.vector()>0.1 does not make much sense in this context, I think you mean to evaluate the function h0 at the coordinate x.

class boundary(SubDomain):
  def __init__(self,h0):
    self.h0=h0
    SubDomain.__init__(self)
  def inside(self,x,on_boundary):   
    return (x[0] < 0.1)and(self.h0(x)>0.1) and on_boundary

---

Thank you very much for your quick and efficient answer.
yes, you totally understood what I meant.
I just had to change the line

bc= DirichletBC(V,max(h0.vector()),BC,'pointwise')

and consequently erase the "on_boundary" from the boundary definition

  def inside(self,x,on_boundary):   
return (x[0] < 0.1)and(self.h0(x)>0.1)

and it worked like a charm !