"Found no facets matching domain for boundary condition" for easy Stokes on Unit Square case

Question

Hi everybody,

I am stuck all the day with the (supposedly) very easy Leaky Cavity - Stokes problem.
I just created create the 2D UnitSquare mesh.

UnitSquareMesh mesh(96, 96);

and, I needed to define the:
a) the lower and side edges where a fluid at rest condition has to be imposed
b) the upper edge where a fluid sliding with a unitary x speed is imposed
c) the upper/right corner, where the pressure is set to zero (being in the stabilized formulation the pressure not defined.

I just write here the classes for the BCs:

class Dirichprex_boundary : public SubDomain
{
    bool inside(const Array<double>& x, bool on_boundary) const
    {
        return  on_boundary && ( near(x[1],1.0) && near(x[0],1.0) );      //  Top/Right Corner x,y = 1
    }
};

class Noslip_boundary : public SubDomain
{
    bool inside(const Array<double>& x, bool on_boundary) const
    {
        return on_boundary && ( near(x[1],0.0) || near(x[0],0.0) || near(x[0],1.0)    ); // left,right and bottom edge
    }
};

class Sliding_boundary : public SubDomain
{
    bool inside(const Array<double>& x, bool on_boundary) const
    {
        return  on_boundary && near(x[1],1.0) ;    //   Top Edge y = 1
    }
};

These are the associated functions:

class Noslip : public Expression
{
public:
    Noslip() : Expression(2) {}
    void eval(Array<double>& values, const Array<double>& x) const
    {
        values[0] = 0.0;
        values[1] = 0.0;
    }
};

class Sliding : public Expression
{
public:
  Sliding() : Expression(2) {}
  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 1.0;
    values[1] = 0.0;
  }
};

In the main, I first define the spaces on the mesh,

Stokes::FunctionSpace W(mesh);
SubSpace Wv(W, 0);
SubSpace Wp(W, 1);

then i define the boundaries and functions, creating thus the DirichletBC.

Noslip_boundary noslip_b;
Sliding_boundary sliding_b;
Dirichprex_boundary dirichprex_b;

Noslip noslip;
Sliding sliding;
Constant zeroprex(0.0);

DirichletBC dirichprexBC(Wp, zeroprex, dirichprex_b);
DirichletBC noslipBC(Wv, noslip, noslip_b);
DirichletBC slidingBC(Wv, sliding, sliding_b);

Finally, I collect the vector with the values:

  std::vector<const DirichletBC*> bcs;
  bcs.push_back(&noslipBC); bcs.push_back(&slidingBC); bcs.push_back(&dirichprexBC);

which I use for defining (and solving) the variational problem:

LinearVariationalProblem problem(a, L, w, bcs);

The results are great in terms of the velocity, whereas the pressure distribution is good, but it is not zero on the wanted vertex (so the pressure is shifted).
In fact, at run-time I get:

Solving linear variational problem.
  *** Warning: Found no facets matching domain for boundary condition.

Do you have any ideas about where the "error" is?
Thank you in advance,
R.
P.S. The code is based on the DEMO C++ Stokes-stabilized, with the difference that I am not importing but creating the mesh and the subdomains.

Answer 1

Your problem is likely due to c), since there are no facets that can possibly match a corner. You need to use a pointwise DirichletBC. For other options on how to deal with the pressure issue, see, e.g., this link.

Comments

Thank you Mikael,
that makes sense. I thought the BC where specified on the vertices. But it make now sense that, if the above methods of declaration look for facets, a corner cannot be recognized.
Do you have any fast suggestion (or example) where I can learn a little bit more about pointwise Dirichlet BC?

---

Just look at the docstring

pydoc dolfin.DirichletBC

and drop the on_boundary parameter.

---

Thank you Mikael.
I was able to solve the problem thanks to your advices.
Just for somebody that may have the same problem in the future, I report below the changes I performed.
On the Subdomain class, I dropped the on_boundary &&

class Dirichprex_boundary : public SubDomain
{
    bool inside(const Array<double>& x, bool on_boundary) const
    {
        return  ( near(x[1],1.0) && near(x[0],1.0) );      //  Top/Right Corner x,y = 1
    }
};

On the main, I just added the pointwise option to DirichletBC.

DirichletBC dirichprexBC(Wp, zeroprex, dirichprex_b ,  "pointwise" );

R.