*** Warning: Found no facets matching domain for boundary condition.

Question

Hello,
I've read the already answered questions on this topic, but I still do not understand why I get this warning. Can anyone help?

Here is an extract of my code. To give you a little background, I'm trying to compute the solution of the time dependent acoustic wave equation on a 2D rectangular domain:

div(grad(u)) - 1/c^2 d^2u/dt^2 = f

I choose f(x,y,t) so that the solution u(x,y,t) = sin(x)/(1+t) -- I am interested in checking the convergence of the numerical solution, knowing the exact solution

I only report here the very first steps (without the time loop).

from dolfin import *
from mshr import *
import numpy as np



size = 20

# Create mesh 
domain = Rectangle(Point(0., 0.), Point(pi, pi))
mesh = generate_mesh(domain, 20)

# Newmark Constants
beta_ = 0.25
gamma_ = 0.5

# Time Constants
dt =  1./size
T = 1.
t = 0.


# Define Function space, trial and test functions
V = FunctionSpace(mesh, 'CG', 1)
a, phi = TrialFunction(V), TestFunction(V)


# Fields from previous time step (displacement, velocity, acceleration)
u0IC = Expression('sin(x[0])', degree=1)
v0IC = Expression('-sin(x[0])', degree=1)

u0 = interpolate(u0IC, V)
v0 = interpolate(v0IC, V)
a0 = Function(V)
a1 = Function(V)

# Boundaries
class Horizontal(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - pi) < DOLFIN_EPS))

class Vertical(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[0] - 0.) < DOLFIN_EPS) or (abs(x[0] - pi) < DOLFIN_EPS))


boundaries = FacetFunction("size_t", mesh, 0 ) # this index 0 is an alternative to the command boundaries.set_all(0)

Horizontal().mark(boundaries, 1) 
Vertical().mark(boundaries, 2) 

bcV = DirichletBC(V,Constant(0.),2)


# Define time-dependent source
c = Constant(0.5) # velocity
source = Expression('( -sin(x[0])/(1.+t) - 2.*sin(x[0])/(pow(c*(1.+t),2)) )', c= c, t=0.0, degree=1)


# Left hand side
m = 1/pow(c,2) * a * phi * dx 

s = inner(grad(a),grad(phi)) * dx - a.dx(0) * phi * ds(2)
s0 = inner(grad(u0),grad(phi)) * dx - u0.dx(0) * phi * ds(2)

f = -source * phi * dx

A = assemble(m + beta_*pow(dt,2)*s)


# Right hand side
b0 = assemble(f-s0)

# Apply BC
bcV.apply(A,b0) 

# Compute a0
solve (A,a0.vector(),b0)

# Save a0 and u0
vtk_file2 = File("ex_sol_results/pressure_acc.pvd")
vtk_file2 << a0

vtk_file2 = File("ex_sol_results/pressure.pvd")
vtk_file2 << u0

Answer 1

Hi, your are creating the DirichletBC object in a wrong way. Consider the following instead

from dolfin import * 
from mshr import *

size = 20

# Create mesh 
domain = Rectangle(Point(0., 0.), Point(pi, pi))
mesh = generate_mesh(domain, 20)


class Horizontal(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[1] - 0.) < DOLFIN_EPS) or (abs(x[1] - pi) < DOLFIN_EPS))

class Vertical(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ((abs(x[0] - 0.) < DOLFIN_EPS) or (abs(x[0] - pi) < DOLFIN_EPS))


boundaries = FacetFunction("size_t", mesh, 0) 

Horizontal().mark(boundaries, 1) 
Vertical().mark(boundaries, 2) 

V = FunctionSpace(mesh, 'CG', 1)
for i in (1, 2):
    count = sum(1 for _ in SubsetIterator(boundaries, i))
    print 'Facet marked as %d is %d' % (i, count)
    # Space, bdry value, facet function, facet marker value
    assert len(DirichletBC(V, Constant(0.), boundaries, i).get_boundary_values()) > 0

Comments

Thank you MiroK!

So to apply homogeneous Dirichlet BC to a portion i of the boundaries one always needs to do

bcv = DirichletBC(V, Constant(0.), boundaries, i)

Defining explicitly ds = Measure("ds", subdomain_data=boundaries) as chris_richardson was saying is then not necessary?

---

Correct, there is no need for defining ds in order to get DirichletBC work properly.

Answer 2

In your form, you need to define what "ds(2)" means, e.g. with:

ds = Measure("ds", subdomain_data=boundaries)

Comments

Hi,

thank you for your reply. Unfortunately adding this line doesn't make the warning disappear. Do you have other ideas??