Solving with interior boundary conditions

Question

Dear all,

I am trying to solve a boundary problem. Thanks to the answer to a previous question of mine, I now know how to apply a boundary condition in the interior of the domain to a vector:

from dolfin import *
mesh = RectangleMesh(0, 0, 1, 1, 10, 10)
V = FunctionSpace(mesh, "CG", 1)
# initialize mesh connectivity so that Facet.exterior() works
mesh.init()
# define the interior of the domain by looking at each facets
facet_domains = FacetFunction('size_t', mesh)
facet_domains.set_all(0)
for f in facets(mesh):
  if any(ff.exterior() for ff in facets(f)):
    facet_domains[f] = 1
bc_in = DirichletBC(V, Constant(-1.0), facet_domains, 0)
u = Function(V)
bc_in.apply(u.vector())

Unfortunately, I am now trying the same approach on a linear variational system defined on the boundary of the domain and it does not work:

from dolfin import *
mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "CG", 1)
# initialize mesh connectivity so that Facet.exterior() works
mesh.init()
# define the interior of the domain by looking at each facets
facet_domains = FacetFunction('size_t', mesh)
facet_domains.set_all(0)
for f in facets(mesh):
    if any(ff.exterior() for ff in facets(f)):
        facet_domains[f] = 1

bc_in = DirichletBC(V, Constant(-1.0), facet_domains, 0)

u, v = TestFunction(V), TrialFunction(V)
a = u*v*ds
L = Constant(1.)*v*ds

w = Function(V)
solve(a == L, w, bcs=bc_in)

plot(w, axes=True)
interactive()

This script produces the following message in the terminal:

Solving linear variational problem.
  UMFPACK problem related to call to numeric
  *** Warning: UMFPACK reports that the matrix being solved is singular.
  UMFPACK problem related to call to solve
  *** Warning: UMFPACK reports that the matrix being solved is singular.

And here is the resulting plot:

It looks like bc_in has not been applied on the whole interior but only on a set of facets immediately below the domain boundary... I guess I am still missing something important in the way boundary conditions and facets are working... Can somebody help me please ? Thanks !

Answer 1

I'm not sure what is the correct, clean solution but I know the cause. The problem is that sparsity pattern of matrix resulting from bilinear form a does not contain any interior DOFs and setting these rows to identity while applying BC fails. (With uBLAS backend this fails silently which is possibly a bug. With PETSc error message is issued.)

For me working workaround is

a = u*v*ds + Constant(0.0)*u*v*dx

(Note that 0.0*u*v*dx would not work as UFL optimizes multiplication by zero out.) However this may not work with all versions of DOLFIN/PETSc as PETSc may optimize assembled zeros out.

Another option is to use keep_diagonal=True assembler instruction.

assembler = FooAssembler(...) # some assembler
# this instructs assembler to keep diagonal entries in sparsity pattern
assembler.keep_diagonal = True 
tensors = assembler.assemble(...)

# or we can use assemble shorthand
A = assemble(..., keep_diagonal=True)
# or
A, b = assemble_system(..., keep_diagonal=True)

# and solve linear algebraic system
solve(A, solution.vector(), b)

keep_diagonal=True should work from 1.2.0 release.

Comments

Again, this is quite easy, I just forgot about DOFs at vertices. See comment above for fix.

---

Ok, now it indeed works. Thank you very much !

Now, out of curiosity (since an order of 2 is fine for me now), what would you suggest if the function space was of order 3 or more ? Then midpoints are not enough !?

---

In principle same procedure. But there are no DOFs at facet midpoints. Just type

dolfin-plot Lagrange triangle 3

to your command-line to see where are appropriate DOFs located.

---

Ok, thanks. One last question (sorry!) : What is the proper method to get the endpoints of the facet ?

---

Try

 f = Facet(mesh, 0)

 for v in vertices(f):
     p = v.point()
     # do something with p

If you're considering programming this for arbitrary order CG space, you may appreciate member functions of V.dofmap().

But the most robust way is probably to avoid DirichletBC at all and manipulate DOFs of FE tensors directly. You can probably avoid "topology -> DOF coordinates -> DOF index" pipeline and skip directly to "topology -> DOF index".